ACULSR.036A 



PATENT 



IN THE UNITED STATES PATENT AND TRADEMARK OFFICE 



Applicant 



Jackson Streeter et al. 



Group Art Unit 3739 



Appl. No. 



10/764,986 



Filed 



January 26, 2004 



For 



LOW LEVEL LIGHT THERAPY FOR 
THE ENHANCEMENT OF 
NEUROLOGICAL FUNCTION 



Examiner 



Henry M. Johnson III 



DECLARATION OF LUIS DE TABOADA PURSUANT TO 37 C.F.R. S 1.132 



Commissioner for Patents 

P.O. Box 1450 

Alexandria, VA 22313-1450 



I, Luis De Taboada, declare as follows; 

1. I am a co- inventor of the claimed subject matter of the above-captioned patent 
application. 

2. I have reviewed the above-captioned patent application, including the 
specification and figures, the claims as originally filed, and the pending claims as amended in the 
"Amendment and Response to May 31, 2007 Office Action" submitted herewith. I have also 
reviewed the May 3 1 , 2007 Office Action in the above-captioned patent application, and the prior 
art reference cited therein, namely, PCT Publication No. WO 99/62599 to Oron ("Oron"). 

3. At page 6, lines 21-27, Oron discloses invasive phototherapy experiments that 
were performed on rats by cutting an incision in the scalp to expose the skull and attaching a 
laser source to the parietal bone in the skull. Oron further discloses that an energy density of 8 
mW per square centimeter reached the brain tissue of the rats. At page 7, lines 7-9, Oron further 
discloses phototherapy experiments that were performed on white mice by irradiating through the 
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skin that covers the skull. Oron further discloses that an energy density of 5 mW per square 
centimeter reached the brain tissue of the mice. 

4. A rat or mouse brain is typically only about 10 to 12 millimeters in diameter, so 
no portion of the rat brain is at a depth of approximately 2 centimeters below the dura of the 
brain. In addition, due at least in part to the drastically different size of a human brain versus a 
rat or mouse brain, it is not obvious that irradiation of the surface of a rat or mouse brain would 
have a corresponding effect if performed on a human brain. Persons skilled in the art of 
phototherapy would understand the disclosure of Oron as not disclosing or suggesting non- 
invasively irradiating a portion of a brain with light having a power density at least about 
0.01 mW per square centimeter at a depth of approximately 2 centimeters below the dura. 

5. At page 8, lines 4-13, Oron also discloses measurements of the penetration of 
laser irradiation through a fresh human skull to illustrate the transmittance of light through the 
skull wall. Oron discloses that the beam diameter of the laser was dispersed from 2 millimeters 
externally to the skull to 3.5 centimeters in the skull cavity and had a power density of 3 mW per 
square centimeter after penetration through a skull wall having a thickness of about 
8 millimeters. Oron is silent regarding the transmission of light through human brain tissue. 
Oron is also silent regarding the irradiation of the brain with light at a point approximately 
2 centimeters below the dura. Persons skilled in the art of phototherapy would understand from 
the disclosure of Oron that Oron does not disclose or suggest irradiating a portion of the brain 
approximately 2 centimeters below the dura with light having a power density at least about 
0.01 mW per square centimeter. 

6. Transmission of light through a material is described by the Beer-Lambert Law, 
which can be expressed as l{L)={l-Rp)*Io*e~^'^ , where lo is the input intensity, i is the 
linear distance from the input point to the measured point, I(L) is the intensity of light after 
propagating through a thickness L of the material, Rf is the relative mean refractive index of the 
material, and fi, is the total attenuation coefficient. The total attenuation coefficient can be 
expressed as: //, where /a^ is the absorption coefficient, is the scattering 
coefficient, and g is the scattering anisotropy parameter. The quantity {\-g)fis is referred to as the 
reduced scattering coefficient fj. V 
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7. The published literature on the optical properties of brain tissue at the time of 
Oron stated that the absorption coefficient ^la of brain tissue was about 0.090 cm"', and the 
reduced scattering coefficient n's of brain tissue was about 8.2 cm"^ {See, e.g., Pogue et al., 
"Comparison of imaging geometries for diffuse optical tomography of tissue," Optics Express, 
Vol. 4, No. 8, 12 April 1999, pages 270-286; Bevilacqua et al., "In vivo local determination of 
tissue optical properties: applications to human brain," Applied Optics, Vol. 38, No. 22, 1 
August 1999, pages 4939-4950; Firbank et al., "A Theoretical Study of the Signal Contribution of 
Regions of the Adult Head to Near-Infrared Spectroscopy Studies of Visual Evoked Responses" 
Neuroimage, Vol. 8, 1998, pages 69-78; each of which is included in the attached Exhibit A.) 
Assuming that the relative mean refi-active index of brain tissue is equal to zero (thereby 
providing the maximum transmitted intensity), the Beer-Lambert Law would be used by persons 
skilled in the art to determine that application of a power density of 8 mW/cm^ at the surface of 
the human brain would yield about 5 x 10''' mW/cm^ at a depth of 2 centimeters. Thus, persons 
skilled in the art would not expect that application of the power densities disclosed by Oron 
would result in a power density of at least about 0.01 mW/cm^ at a depth of 2 centimeters below 
the dura. 

8. In addition, using these accepted prior art values for the absorption and 
transmission coefficients for brain tissue, persons skilled in the art would expect that in order to 
get a power density of at least about 0.01 mW/cm^ at a depth of 2 centimeters below the dura, a 
power density of about 159 W/cm^ would be needed at the surface of the brain. Persons skilled 
in the art would recognize that this power density would cause significant damage of the patient's 
brain and would not be practical for treatment. 

9. Despite the teachings of the prior art that the transmission of light through brain 
tissue was too low to make phototherapy of the brain practical, the above-captioned patent 
application describes a technique of enhancing neurological function in a patient by 
noninvasively irradiating the patient's brain. To provide support for our contention that the 
actual transmission of light through brain tissue is higher than was expected trom the prior art 
references, I took part in a pair of separate studies in which the transmission of radiation through 
a cadaver human skull and through a cadaver human brain was measured. In the skull 
transmission study (V.V. Lychagov et al, ''Experimental study of NIR transmittance of the 
human skull," Complex Dynamics and Fluctuations in Biomedical Photonics III, ed. by V.V. 
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Tuchin, Proc. of SPIE VoL 6085, 6085'r (2006), attached herewith as Exhibit B), an input beam 
having a wavelength of 81 0 nm, a power density of about 68 mW per square centimeter, and a 
beam diameter of about 30 millimelers was applied to the shaved skin of the skulls of 20 
cadavers at five regions with various thicknesses. Figure 2 of Exhibit B illustrates that the 
transmittance of the light varied fiom 0.5% to 5%, depending on flie thiclcncss of the irradiated 
portion of the skulL We perfonBcd an independent measurement (unpublished) of the 
transmission of 80S-nm Hglit with a power density of 424 roW pear square centimeter applied to a 
cadaver human skull (includirtg tlie shaved skin) having an average thickness of about 
11 millimeters, and calculated the overall transmittance of the skull tissues to be about 2%, 
which agrees well with, the results of Exhibit B, 

10, In the brain transmission study (unpublished), measurements showed that for an 
input beam having a wavelength of 808 nm, a power density of 10 mW per square centimeter 
applied directly to the dura, and a beam diameter of 30 millimeters at the dura, the power 
density at a depth of 1 .5 centimeters irom the dura was 0. 148 mW per square centimeter, and the 
power density at a depth of 2,5 oentimcftors fl-om the dura was 0.01 0 mW per square centimeter. 

11. T liereby declare that all statements made herein of my own knowledge are true, 
and that all statements made upon information and belief are believed to l« Ixue; and flirther, that 
these statements were made with the knowledge that willlul, false statements and the like so 
made arc punishable by fine or imprisonment or both under Section 1001, Title 18 of the United 
States Code, and that willful, false statements may jeopardize the validity of the ^plicatian or 
anypatcut issuing thereon. 

Dated: /oXk"' y^QOy Bv: / -y^:^^.^'^ /C^^tsC^ 
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EXHIBIT A 



Comparison of imaging geometries for 
diffuse optical tomography of tissue 

Brian W. Pogue, Troy O. McBride, Ulf L. Osterberg and Keith D. Paulsen 

Thayer School of Engineering, Dartmouth College. Hanover NH 
pogue@dartmouth. edu 

Abstract: Images produced in six different geometries with diffuse optical 
tomography simulations of tissue have been compared using a finite 
element-based algorithm with iterative refinement provided by the Newton- 
Raphson approach. The source-detector arrangements studied include (i) 
fan-beam tomography, (ii) full reflectance and transmittance tomography, as 
well as (iii) sub-surface imaging, where each of these three were examined 
in a circular and a flat slab geometry. The algorithm can provide 
quantitatively accurate results for all of the tomographic geometries 
investigated under certain circumstances. For example, quantitatively 
accurate results occur with sub-surface imaging only when the object to be 
imaged is fully contained within the diffuse projections. In general the 
diffuse projections must sample all regions around the target to be 
characterized in order for the algorithm to recover quantitatively accurate 
results. Not only is it important to sample the whole space, but maximal 
angular sampling is required for optimal image reconstruction. Geometries 
which do not maximize the possible sampling angles cause more noise 
artifact in the reconstructed images. Preliminary simulations using a mesh 
of the human brain confirm that optimal images are produced from 
circularly symmetric source-detector distributions, but that quantitatively 
accurate images can be reconstructed even with a sub-surface imaging, 
although spatial resolution is modest, 
©1999 Optical Society of America 

OCIS code: (1 1 0.6960) Tomography, {1 70.3660) Light propagation in tissues, 
(170,3010) Image reconstmctiw teclmigues, (iOO.3190) Inverse Problems 
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1. Introduction 

DifRise Optica] Tomography (DOT) imaging of tissue is emerging as a viable method for 
non-invasively determining hemoglobin concentration, hemoglobin oxygen saturation, 
cytochromes, lipids and water in vivo on a centimeter spatial scale. In the past decade, 
promising applications have emerged in areas such as neonatal and adult cerebral 
monitoring [1, 2] and breast cancer diagnosis [3]. In addition, DOT may be even more 
promising when combined with traditional imaging methods [4-7] or when coupled to 
exogenous contrast agents to monitor tissue fiinction, either through fluorescence or 
absorption based imaging [2, 8, 9]. One of the limitations in this field of research has been 
that the rapid advances in technological instrumentation have eclipsed the development of 
accurate image reconstruction algorithms for diffuse tomography. In this study we focus on 
one particularly useful image reconstruction algorithm, and apply it in several different 
imaging configurations in order to evaluate how the source-detector- tissue geometry affects 
the resulting image. 

The use of an appropriate reconstruction algorithm is important because the inverse 
problem for DOT is inherently ill-posed, and in many cases is ill-conditioned [10], so that 
standard linear approaches provide limited resolution and contrast inaccuracy [11], Several 
iterative reconstruction methods have been tested, and in particular, the Newton-Raphson 
approach has been shown to provide quantitatively accurate images of the absorption 
coefficient in a circular tomographic geometry [12-16]. This approach has been developed to 
characterize breast tumor optical properties in order to diagnose changes in hemoglobin 
concentration and oxygen saturation of normal and malignant tissues, and is currently being 
used in our clinical prototype imager. 




Fig. 1. Schematic of three typical source-detector geometries used in diffuse 
optical tomography applications. The projection-shadow geometry (a) is used in 
several optical mammography scanners [3, 23-25], while the circular-tomography 

geometry (b) is used for both breast and brain imaging [16, 26], as is the sub- 
surface imaging geometry (c) [2, 27, 28]. 



There has been considerable study of the type of measurements needed to provide 
optimal sensitivity in separating the scattering and absorption coefficients in tissue [17-19], 
There has also been considerable study of the resolution limits of diffuse tomography [20- 
22]; however, in general the optimal data types and the resoluton limits are liinctions of the 
reconstruction algorithm used to retrieve the images and the source-detector arrangement. 
Resolution limits calculated from forward solutions may not apply to inversion algorithms 
which preserve the non-linearity between the internal optical properties of interest and the 
measured projections at the boundaries. Similarly there has been little investigation of the 
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geometry and how its boundaries affect the imaging process. As a practical matter, it is quite 
important to assess the effects that boundaries and source-detector configurations have on the 
image reconstntction in order to design prototype medical systems. These design 
considerations are ideally assessed through computer simulations to avoid costly and wasteftil 
construction of clinical prototype systems which may have inherent design flaws. With this 
goal in mind, we assess geometries which to date represent the most common configurations 
used for near-infrared imaging and spectroscopy. 



nj ilAi-: tClX-C IVlK'f! 



Fig. 2. Geometries used in this study are shown, where the outward arrows indicate 
(ielectoT luciilums and the inward arrows denote source locations. The geometries 
include (a) circular reflectance and transmittance imaging which uses circular 
symmelry for sampling angles (b) circular fan-beam imaging (where the source- 
detector array is rotated around the tissue), (c) sub-surface tomography which only 
uses a sub-set of projections irom the surface of the tissue, (d) flat slab reflectance 

and transmittance imaging using maximal sampling from both surfaces, (e) flat 
slab fan-beam which does not use any reflectance data, and (e) sub-surface imaging 
with a flat upper surface. Together these six geometries, in various forms, define 
the possible configurations for imaging tissue with DOT. 
The geometries chosen for many laboratory and clinical near-infi-ared imaging studies 
fall into three general catagories, which are delineated here as (a) projection-shadow, (b) 
circular-tomography and (c) subsurface imaging, as shown in Fig. 1. Each of these 
geometries is being used in clinical investigations for imaging either breast or brain tissues 
with near-infrared light; yet, there has been little comparison of how the source-detector 
array geometry will alTects tlie quality of the optical measurements and resulting image. It 
should be noted that the projection-shadow geometry is generally used for traditional x-ray 
film style images, where there is no inversion of the data to compute tomographic images 
[3], and that this geometry does not lend itself to tomographic reconstruction. However, even 
with this restriction, Fantini e( al. have been able to retrieve quantitative optical properties 
firom breast tumors using the assumption that the tumor is roughly spherically shaped [29], 
and this procedure can be even more accurate when multiple projections and multiple 
wavelengths are incorporated into the analysis [24]. Similarly, Quaresima et al. have used 
this geometry to measure the spectral signature of breast tissue [30], A variant of this 
geometry, which exploits a single source location and the angular projections to multiple 
detectors, has been employed by Li et al. [31] to produce backprojection images of breast 
tissue. Circular tomography has been used successfully in our own studies of hemoglobin 
based imaging [15, 16] which is being developed for breast cancer diagnosis [32], as well as 
in work of other investigators for cerebral imaging [26, 33]. 
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In this study, the same DOT image reconstruction algorithm is used to examine: (a) 
circular full reflectance and transmittance tomography (b) circular fan-beam tomography, (c) 
sub-surface imaging under a curved boundary, (d) flat slab reflectance and transmittance 
tomography, (e) flat slab M-beam tomography, and (f) sub-surface imaging below a flat 
surface, as shown in Fig. 2. The goal of these tests is to examine the influences of (i) the 
tissue geometry and (ii) the source-detector configuration upon the recovered image accuracy 
and quality, and to determine the optimal data collection strategy for characterizing tissue 
optical properties with this algorithm. The fiat slab geometry (Fig. 2 (d) to (f)) is well suited 
for breast imaging since it provides a configuration similar to x-ray mammography, but it is 
not clear how the image quality will compare to the more symmetric sampling provided by 
the circular geometries (Fig. 2 (a) to (c)). We also compare a fan-beam orientation of the 
detectors where measurements are only recorded on the opposite side to the source (Fig. 2 (b) 
& (e)), with the full set of projections from all points equally spaced around the tissue (Fig. 2 
(a) & (d)). The former orientation is technologically easier to implement because alt 
measurement intensities are within the same range of detector sensitivity, so that parallel 
implementation is more easily facilitated. However, it is not obvious how this decrease in 
data affects the resulting optical property distribution, as compared to the full reflectance and 
transmittance sampling. Finally, subsurface tomography is examined since this is a geometry 
which is being widely used for near-infrared tomography of brain and breast tissues [2, 28, 
34]. 

2. Theory & Computational Methods 

The theoretical details of diffuse optical tomography have been described by several authors, 
and some good overviews are included in references [35-37]. Forward calculations of diffuse 
light in tissue can be used to simulate experimental measurements of heterogeneous tissue, 
and by matching calculated measurements to these simulated experimental measurements an 
image of the optical properties can be recovered. Since frequency-domain light 
measurements are used in this work, the frequency-domain difiiision equation is solved with 
a finite element method. The general form of this equation is: 

V ■ D(r)V<D(r, £0) - [^(r) + 1(0) I cj\^{r, (o) = S<,{r, a) 

(1) 

which describes optical fluence rate ^{r,m) at position r, and frequency to, given the 
distribution of absorption coefficient, ytjj) and diffusion coefficient D(r), where c is the 
speed of light in the medium^ and So(r,co) is tlje optical source. The diffusion coefficient is 
usually defined as D = [3(jis+l^a)]''. where |Js is the transport scattering coefficient of the 
medium. This equation is solved with Type Hi (Robin) boundaiy condition on a finite 
element mesh, where all exterior elements are subject to: 

O(^) + AV(t(^,£(j)-o = 0 

where ^ is the position on the boundary, h is the unit vector normal to the boundary, and k is 
a calibration constant which has been determined empirically by matching the forward 
calculations to well controlled experiments [35], From this solution, we calculate a set of 
optical measurements, which in this case are phase shift at frequency, o), and the 
logarithm of the signal amplitude. While the solution to equation (1) is straightforward 
given the distribution of and n/, recovering the image distribution of fx, and \lI requires 
an iterative solution, since equation (1) cannot be directly solved for these values, in part 
because 0(r,c£)) is not known inside the tissue, but only at the discrete points (detector 
locations) on the boundary, 

The Newon-Raphson iterative method for finding the solution of Hs and \4 starts by 
solving for the minimum of the functional 
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z(n)=i<t>'-<^''r(<i>'-4>°)+^(ti) 

(3) 

where we use |a to represent both ^a and p./ as the set of coefficients, and F(n) is an arbitrary 
fiinctional which smooths the solution to reduce the ill-posed nature of the eventual matrix 
inversion [38], Minimizing (3) by taking the derivative and setting it equal to zero, followed 
by expansion of the solution in a Taylor's series [39], leads to, 

(4) 

where the elements of J are given by d^;/d\Xk, for observation i, at each pixel position k. In 
this derivation, the higher order derivatives of (j) are assumed to be negligible, and F(^) is 
chosen to result in the classical Tikhonov regularized formulation [38]. Finally the iterative 
solution becomes, 

(5) 

where is the update for the current estimate of \i. The values of if^ and J are calculated 
using the current estimate of [i. In practice we have found that the regutarization parameter 
can be varied between iterations for improved reconstruction, and can be varied spatially to 
compensate for some geometrical artifacts induced during the reconstruction process [15]. 

In the following simulations, the forward and inverse calculations have been carried out 
with a finite element solution which has been tested extensively in previous papers [ The 
simulated measurements were calculated on the same meshes that were used for the 
reconstructions and 1% random gaussian-distributed noise was added to both the real and 
imaginary components of the measurements in order to simulate true source-detector noise. 
All reconstructions were carried out with the same regularization parameter, and 5 iterations 
have been used in each case. Each iteration requires approximately 1.5 minutes on a 200 
MHz desktop computer running Windows NT. In all the reconstructed images shown in this 
paper, the dimensions of the field are shown in millimeters on the x and y axes of the figure. 

3. Results 

3. 1 Imaging Fields with a Single Localized Target 

A series of forward and inverse simulations have been completed for each geometry shown in 
Figure 2. The overall goal has been to maintain a uniform test between geometries and 
determine those geometries which introduce significant artifacts into the reconstruction by 
the nature of the boundaries and source-detector configuration. There are an infinite number 
of possible tests which could be used in each of these geometries, but we have chosen to 
begin with a physiologically relevant contrast in a tumor sized target. In all cases the 
background optical properties were fixed at Ha= 0.01 mm"' and \xl = 1.0 mm"', and a 
sunulated target within the field has been included with a radius of 10 mm and optical 
properties of = 0.02 mm"' and n/ = 1.0 mm"'. In each geometry the target is moved 
around to positions within the imaging field which have different levels of symmetry, 
potentially affecting the resulting image qualtity. 

3. 1. 1 Circular Reflectance and Transmittance Imaging 

A circular mesh was used with 221 nodes in the coarse calculation of the Jacobian and 
matrix inversion as shown in equation (5), and 425 nodes in the fine mesh were used for 
solution of equation (1). A circularly symmetric configuration of 16 sources and 16 detectors 
was used in the solution, providing a data set of 256 projections, as shown in Fig. 2 (a). The 
field was 86 mm in diameter, with the sources and detectors spaced symmetrically around 
the periphery, with one every 11, 25 degrees alternating between source and detector. The test 
target was moved to different radial positions between the center and the edge of the field to 
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examine the quality of the reconstructions which are presented in Fig. 3. Because this 
simulation has circular symmetry, no other target locations are needed to characterize the 
image field. ______________ 




Fig. 3. Reconstructed set of images (video sequence) of the test target at different 
vertical positions within the circular mesh, using suurcts-detector geometry in Fig, 2 

(a) with a circularly symmetric distribution of 16 source and 16 detectors for 
reflectance and transmittance imaging. The color bar on the right is in absorption 

coefficient units of mm"'. 

3.1.2 Circular Fan-Beam Imaging 

This orientation of source and detectors is shown in Fig 2(b), where for each source location, 
8 detectors were used on the opposite side of the phantom. Then the source and locations of 
the 8 detectors are rotated around the phantom to 16 different sets of projections through the 
phantom, resulting in 128 measurements of phase shift and signal amplitude. The spacing of 
field and source-detectors was otherwise identical to Fig. 3. This geometry has the 
experimental advantage that the signal intensity at all 8 detectors is within an order of 
magnitude in variation, so that large changes in the dynamic range of the detectors is not 
encountered. This data set can be considered as a sub-set of the data available in the 
previous case of fliil reflectance and transmittance with the near-source detectors removed. 
The location of the test target was varied in the object field between the center and the top 
edge. Again, since the sources and detectors are located in circularly symmetric orientations, 
this set of target locations is sufficient to characterize the entire image field for this size and 
shape of object, as shown in Fig. 4. 




Fig. 4. Reconstructed set of images (video sequence) showing the test target at 
different vertical positions within the circular mesh, using the fan beam geometry 
with 16 sources transmitting to 8 detectors. The color bar on the right is in 
absorption coefficient units of mm'', 

3.1.3 Sub-surface Imaging- Curved Boundary 

Near-infrared imaging is confounded by the high attenuation of the light in tissue, such that 
in practice, it is often not possible to detect light reliably through tissue samples thicker than 
10 cm. Also, the instrumentation for imaging through tissue less tiian 5-6 cm can be much 
simpler, more robust and less expensive than for thicker tissues [28]. For this reason, many 
investigators have focused on sub-surface imaging applications where the sources and 
detectors are located on the same tissue surface and are arranged to probe the interior of the 
tissue. This sub-surface tomography can be performed on a curved tissue such as a neonate 
or adult head and female breast. These tissues are well simulated with a finite element 
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solution assuming that the boundary curvature is known and can be input into a finite 
element mesh. If the light does not sample deeply into the tissue, the tissue boundary can be 
approximated by a simple circular mesh which has a matching radius of curvature. In the 
next series of simulations, this sub-surface imaging method has been investigated using the 
same circular mesh as above, but with only the source and detector locations on the upper 
region of the mesh activated. Both situations illusfraled in Fig. 2 (c) were examined 
includmg: 8 sources-8 detectors (8x8) and 4 sources-4 detectors (,4x4), The imaging field is 
symmetric about the center of the source-defeclor array, so that the image quality is 
adequately assessed by translating the test object vertically down from the upper surface, as 
shown in Fig. 5. 



Fig, 5. Reconstructed setof im^es of the tusl tarset within a circular mesh using 
Ihe sub-surface orientation with (a) 8 souitos imd S dLicuioi s m an arc, and (b) 4 
sources and 4 detectors in a smaller arc on llie siii luce (as slmw n in hig. 2(c)). 
Color bar at nght ism absorption cocllicicni unil.s oi mm '. 

3. J. 4 Slab Reflectance and Transmittance Imaging 

The slab geometry mimics x-ray mammography and minimizing the thickness of tissue that 
must be imaged, potentially allowing the use of less sensitive, but more stable, solid state 
detectors. Light signals from the tissue can be recorded from both the upper surface near the 
source, as well as transmitted through the slab. This geometry is referred to here as slab 
reflectance and transmittance imaging, and is examined in the same set of simulations as the 
circular geometry described in the previous sections. One major difference from the circular 
geometry is the lack of symmetry, so that the response of the imaging field will likely vary 
with both vertical and horizontal position. Thus, in order to fully characterize the image 
plane, the test object has been translated along three lines: vertically through the center, 
horizontally through the center, as well as horizontally along the upper surface, as shown in 
Fig. 6, In the slab mesh, 260 nodes were used in the coarse calculation of the Jacobian in 
equation (5) and 380 nodes were used in the line mesh solution of equation (1), The source 
detector arrangement was as shown in Fig. 2(d), with 5 mm between each source and 
detector. 
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Fig. 6. Reconstructed set of images of the test target at different (a) vertical 
positions (b) lateral positions along the central line, and (c) lateral positions along 
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the surface, using the slab reflectance and transmittance geometry of sources and 
detectors as shown in Fig. 2 (d). 

3.1.5 Slab Fan-Beam Imaging 

The fan-beam geometry for slab breast imaging is an obvious choice since it can involve the 
use of parallel detection through a CCD or other parallel detectors, which are highly stable 
and provide robust detection. The simplicity of translating the source across the upper 
surface of the tissue is beneficial, but it is not clear how the loss of reflectance information 
will affect the reconstructed images of the tissue interior. The geometry used in this section 
was as shown in Fig 2.(e), with 5 mm between neighbouring fibers. In these simulations, the 
asymmetry is simitar to the previous section, so that the same set of simulations were carried 
out to examine the imaging field sensitivity, shown in Fig. 7. 




Fig. 7. Reuoiislructed set of images ol the te.st target at different (a) vertical 
positions (b) lateral positions along the central line, and (c) lateral positions along 
the surface similar to those in Fig. 6, using tan-beam slab imaging. Color bar at 
right is in absorption coefficient units of mm . 



3. J. 6 Sub-surface Imaging- Flat boundary 

The flat boundary is perhaps one of the most often used geometries for near-infrared sub- 
surface imaging because it is simple, readily available clinically in many body locations, and 
lends itself to analytic reconstruction methods. This geometry was tested here with 16 
sources and 16 detectors all located on the upper surface of the tissue volume producing 256 
projections through the medium, spaced alternating at in Fig. 2 (f), with 2,5 mm between 
each. The high number of sources and detectors was employed to provide the best 
performance that we could obtain in this geometry. The test object was translated vertically 
down from the surface, as with the sub-surface imaging on the circular mesh (Section 3. J. 3), 
with the results shown in Fig. 8. 




Fig. 8. Reconstructed set of images where the test target was translated verticafly 
from the middle of the tissue volume (25 mm down) up to the surface, using 16 
sources and 16 detectors alternated and equally spaced along the upper surface of 
the region, as illustrated in Fig. 2 (f). 

3.2 Comparison of Geometries ~ Contrast and Noise 
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The above simulations do not completely define the characteristics of the imaging geometry 
and slight improvements can be achieved for individual cases by adjusting the regularization 
parameter both during an iteration and spatially-distributed across the imaging plane [15]. 
Nonetheless comparing the reconstructed peak values in the target zone and the typical noise 
in the homogeneous background of the image are useful measures of the ability to reconstruct 
objects accurately. In Fig. 9 (left-top graph) the reconstructed peak values for each target are 
plotted for changes in vertical position through the image plane. Fig. 9 (right-top graph) is a 
plot of the standard deviation of the noise within what should be the homogeneous regions of 
the resulting images. In the peak absorption coefficient, if the value is within 20% of the 
original target it is considered acceptable. This criterion is simply based upon the typical 
variations in reconstructed properties observed in previous studies [16]. 




Target daplh from uppor turtaca [mm] Target dapth from upper suitar:B [mm] 



Fig. 9. Calculated (left-top graph) peak absorption coefficient values for the object 

at different positions within the simulated phantoms, and (right-top graph) 
reconstructed noise in the background region of the image measured by sampling of 
random regions of interest outside the target zone. The U'ue value in the first graph 
is shown with a solid horizontal line, and the 20% region of acceptability is defined 
by the dotted lines. Calculated (left-bottom graph) target location as a depth from 
the upper surface, and (right-bottom graph) the calculated target full width at half 
maximum (FWHM) height as a ftinction of depth in the medium, here taken as an 
average of the lateral and vertical directions. 

A comparison of the three subsurface imaging geometries is shown in Fig. 10, using the 
same criteria of acceptable peak reconstruction value, as shown in Fig. 9. The reconstructed 
peak absorption coefficient values for the single targets are shown in Fig. 10 (left-top) as a 
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function of the depth into the medium, and typical noise values are plotted in Fig. 10 (right- 
top). In this case, the only geometry which was able to reconstruct accurate absorption 
coefficient values over a moderate range of depths was with the 8x8 circular arc and the 
reconstructed peak values decrease significantly with depth for all tiiree sub-surface imaging 
geometries. 





Targst dsptti fram upper surface Imm] Target depBl tiom upper surface |ram | 

Fig, 10, Calculated (left-top graph) peak absorption coefficient values for the target 
zone at different positions within the simulated phantoms, and (right-top graph) 
reconstructed noise in the homogeneous background region of the image outside 
the target zones. Calculated (left-bottom graph) target location as a depth from the 
upper surface, and (right-bottom graph) the calculated target full width at half 
maximum (FWHM) height as a function of depth in the medium, measured here 
both in the lateral and vertical directions. 

3.3 Imaging Fields with Multiple Targets 

In addition to imaging single objects it is important to recover more complicated 
distributions, so an additional test has been constructed to compare multiple targets in the 
circular and slab geometries. The test field contained three targets positioned at different 
depths, as illustrated in the circular and slab schematic images in the left column of Fig, 11, 
In this case, two targets were 20 mm in diameter with one in the exact center of the field and 
one on the right edge at the surface (and upper region in the slab case). The third target had 
a 10 mm diameter and was located on the lefil edge of each image (and towards the bottom 
surface in the slab). All three targets presented a contrast of 2:1 in absorption. The resulting 
reconstructions for reflectance and transmittance geometries are shown in the middle column 
while fan-beam reconstructions are contained in the right column of Fig. II, The image 
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quality of the two sub-surface geometries is clearly limited in depth resolution, so that the 
same tests have not been applied to these cases. Therefore, a separate multi-target test 
consisting of two sub-surface objects, both 10 mm diameter and 2:1 in contrast but with one 
on the surface and the other 5 mm below, was used to compare these sub-surface 
configurations, as shown on the lell of Fig. 12. The reconstructed images using the 4x4 and 
8x8 circular arcs and the 1 6x 1 6 flat slab are presented in the right columns of the figure. 




Fig. 1 1 . Rewmslrucleil images of three localized targets within a homogeneous 
background field, (schematics shown in the left column) using geometries of (upper 

middle) circular reflectance & transmittance (upper right) circular fan-beam, 
(lower middle) slab reflectance & transmittance, and (lower right) slab fan-beam. 
The test targets (left) were all 2:1 absroption contrast from the background, and the 
two on the right were 20 mm diameter, while the one on the left is 10 mm in 
diameter. 




Fig. 12. Reconstructed images of the test fields (shown in the left coluni) are 
displayed at right. The test fields contained two 10 mm diameter objects with one 
on the surface and the other 5 mm below the surface. In the top row, the middle 
image shows the result of the 4 source-4 detector array while the right image shows 
the result firom the 8 source-8 detector array. In the bottom row, die result from the 
slab geometry is shown using the 16 source-! 6 detector array. 

3. J Imaging Fields with Layered Structure 

A third test was performed to examine the influence of geometry on the ability to recover a 
layered structure within the medium. In this situation, the field consisted of three concentric 
layers of differing ab.sorjition coefficient with a single circular inhomogenetty located within 
the layers- The schematic of the test object is shown in Fig, 13 in the left column, with the 
circular test field in the upper row, and the corresponding slab test field in the lower row. 
The reconstructed images using the full reflectance & transmittance and fan-beam 
approaches are shown in the center and right columns of the figure, respectively. 
Interestingly, in all cases the circular target could be localized and reconstructed to within 
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20% of the true value, but the layered structure near the surface was best recovered in the 
circular geometry. In particular, the full reflectance and transmittance array provided a 
better reconstruction of the layered structure than did the fan-beam approach, while both slab 
approaches were dominated by artifacts in the background. 

In Fig. 14, the same layered test fields as in Fig. 13 were used to evaluate the sub-surface 
imaging geomefries. These methods were not expected to reconstruct the lower regions of 
the phantom, but the upper layers and the presence of the object can be partially 
reconstructed in all cases. The 4x4 detector strategy and the slab 16x16 scheme cannot fully 
discriminate the object from the layers, while the 8x8 arc geometry can partially reconstruct 




Fig. 13 Tomographic rcconstniclions of a layered tested field with a single target 
included tar both circular (top row) and slab (bottom row) geometries, with two 5 
mm thick layers of Ma = 0.01 and jia = 0,02, with the interior at Ma = 0.015 mm"', 

and a 20 mm diameter object at jia = 0.02 mm~'. The im^es in the left column 
represent the test field, while the middle row of images show the fiill reflectance 

and transmittance images, and the right column of images show the fan-beam 




Fig. 14. Reconstructions of the test fields shown on m the Iclt uulumn using tlie 
sub-surface imaging geometires, including the curved boundary with (upper middle 
image) 4 sources and 4 detectors, (upper right image) 8 sources and 8 detectors in 
an arc, and the slab boundary (lower right image) 16 sources and 16 detectors in a 
flat plane along the surface. 



3.5 Application to Irregular Tissue Geometries - Cerebral Hemoglobin Imaging 
In an effort to consider these source-detector arrangements In a more realistic geometry, a 
magnetic resonance image of a human brain was segmented into a finite element mesh 
(shown in Fig. 15) and used to simulate imaging with the three source-detector geometries (i) 
fan beam (ii) reflectance and transmittance and (lit) sub-surface imaging. The head size was 
scaled to that of a neonate, and typical optical properties were applied to the skin, bone, grey 
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matter and white matter [7,33]. A 15 mm diameter inclusion was added in tlie upper left of 
the brain with Hs = 0.03 mm"', representing a 50% increase from the grey matter, and a 
200% increase from the white matter. The initial guess for the image was a homogeneous 
head shape with [i^ = 0.01 mm"'. The inclusion was observed in all reconstructed images 
with a peak value within 10 % of the true value, indicating that it can be successfully 
characterized with all three source-detector geometries using this algorithm. Interestingly, 
none of the imaging geometries could discriminate the white matter structure from the grey 
matter, but liie fiill reflectance and transmittance method was able to recover the skull optical 
properties accurately, as was the sub-surface imaging geometry. Also the fan-beam approach 
produced significant artifact into the image, as a ringing structure. Surprizingly the fan- 
beam approach recovers the presence of the object clearly, albeit not at a accurate contrast 
level, 




Fig. 15. Reconstructed simulations of a human cranium with the test field (upper 
left), showing absorption coefficients of skin, jia = 0,01, skull, = 0.025, grey 
matter, Ha = 0.02, and white matter, = 0.01 mm"', with a fixed scattering 
coefficient of jj/ = 1.0 mm"'. An artificial inclusion having (*» = 0.03 was included 

in the upper left of the test field. The reconstructed images are shown using fiill 
reflectance and transmittance imaging (upper right), fan beam imping (lower left), 
and 8x8 sub-surface imaging with the detectors in an arc around the upper half of 
the head (lower right). Note diat the sub-surface geometry is not expected to 
recover properties in the lower half of the im^e. 

4. Discussion 

4. 1 Geometry Effecis in the Reconstructed Image 

Determining the optimal source-detector and tissue geometry is not as simple of a task as it 
might first appear because of the interplay between theoretical performance and practical 
constraints, so that some theoretically optimal designs must be considered in light of 
experimental limitations. As a result, all of these issues cannot be fiilly characterized by the 
simple simulations contained in this paper. Also, the reconstruction algorithm could be better 
tailored to any of the specific geometries than in its present generalized form, by altering the 
regularization scheme to improve the reconstruction to provide a more homogeneous solution 
across the imaging field [15]. In the cases shown here, the regularization parameter was 
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chosen to provide optimal reconstructions, through a series of trials, so that the images 
shown here are our best versions. However, even with these uncertainties about 
regularization, there are some reconstruction results which are indicative of the measurement 
data quality which can be used to evaluate and compare the individual geometries. Our goal 
has been to determine if the geometry significantly inhibits the ability to accurately 
reconstruct optical properties, or if the geometry adds artifacts to the resulting image which 
are not easily avoided. 

When comparing reconstruction geometries, it would be ideal to have the same number 
of measurements in the data set, so that changes in geometry are not confounded by changes 
in the conditioning of the matrix; however, in practice the physical size of the source and 
detector fibers can limit the number of possiblilities which are available for a given 
geometry, For example, in the circular geometry, the minimum tissue diameter and the 
diameter of the fibers dictates the number of sources and detectors that are possible. For 
these reasons, we have chosen to compare geometries and source-detector arrangements 
which are physically realistic, while trying to maintain the maximum number of 
measurements in each data set. The fan-beam and sub-surface imaging schemes can be 
thought of as sub-sets of the data measured in the full reflectance and transmittance mode, so 
that it is not surprising that there is some loss of image resolution with each of these two 
methods. We have also constrained the problem to be 2-dimensional in these cases, however 
based upon previous experience we anticipate that similar results are obtained in 3- 
dimensions. Also fully 3-dimensional reconstruction is inherently difficult to interpret, and 
often requires separation into 2-d slices in order to fiilly evaluate the image quality. For 
these reasons, the 2-dimensional images shown here are thought to be representative of the 
available image quality for each respective geometry. 

4.2 Single Targets 

In the video sequences of Figs. 3 to 8, the object is localized accurately in all cases, and is 
close to the original true value at most locations (as shown in Fig. 9). Small deviations from 
the true value can be compensated with changes in the regularization parameters [15], but 
large deviations from the true value which are accompanied by significant noise in the image 
cannot generally be eliminated. An example of this phenomenon is observed in Fig. 7, 
where the fan-beam slab reconstructions show a marked increase in noise within the image, 
while there is no decrease in the recovered peak height of the target. This noise is likely due 
to the decreased angles of projection which are used to sample the field, producing regions 
within the image which cannot be well charaterized by the solution to the inverse problem. 
Maximizing the number of different angular projections through the field is likely the best 
way to improve the accuracy of the inverse solution. 

Recovering images in the sub-surface imaging geometry is more challenging than with 
full tomographic methods, and this is expected since the tissue cannot be sampled as 
effectively ina purely reflectance mode. The most important observation from the data in 
Fig. 9 is that the true value of the object is not easily recovered without being able to sample 
measurements in an arc around the tissue sample, as is the case with the 8 source-8 detector 
array in the circular geometry. By adjusting the regularization parameters, it would be 
possible to reconstruct the peak values for all three geometries used here, as long as the depth 
of the target is not lower than 20-30 mm. At depths below 30 mm, the photon paths cannot 
sample the target effectively. In general as the target is lowered into the medium, the 
reconstructed peak value will decrease. 

4.3 Multiple Targets 

When multiple targets are present as in Figs. 11 & 12, the circular geometries tend to over- 
reconstruct targets on the periphery. This problem can be compensated with a radial 
variation in the regularization parameter [15]. The slab geometries do not recover the 
smaller test target as well as their corresponding circular geometries. There are no clear 
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differences between the fan-beam style versus the full sampling of reflectance and 
transmittance, other than the expected increase in background noise in the fan-beam case. 
The subsurface imaging method suffers firom a field sensitivity which decreases with depth 
into the tissue, such that lower lying objects are not as strongly reconstructed as objects 
immediately at the surface. 

4.4 Layered Fields 

When layered test fields are used, there is serious degradation of the image quality in the 
case of the slab geometry (Fig. 13), whereas both circular geometries were able to provide 
reasonable images of the subsurface ring and the upper inclusion. Again, the fiill reflectance 
and transmittance circular mode provided a slightly more accurate reconstruction of the test 
field than the fan-beam circular mode. Note that the layer in the circular fan beam 
reconstruction is closer to the center than in the original test field (Fig. 13). All three sub- 
surface imaging method are able to reconstruct parts of the layered structure, but only the 
8x8 detection method is sufficient to recover the inclusion quantitatively below the layer. It is 
likely important to ensure that the diffuse projections sample below the suspected inclusion 
in order to produce an accurate reconstruction of the target. This layered test field is 
important, since almost all in vivo tissues have a complex layered structure near the surface 
that may not be clearly resolved unless a circular reflectance and transmittance geometry is 
used. 

4. 5 Cranial Imaging 

For the cranial images in Fig. 15, the complex slrucfure between white and grey matter is not 
well resolved, which agrees with the observations of other investigators [33], The location of 
the skull is well resolved with the full reflectance and transmittance method, but the interior 
is not. The sub-surface imaging geometry provides a method sufTicient to quantify the 
optical properties of the inclusion. The fan-beam imaging appears to image asymmetric 
regions, such as the inclusion quite well but cannot recover the skull or brain details at all. In 
general, the white/grey matter variations which appear on a sub centimeter scale cannot be 
resolved; however, larger changes in the optical properties, such as the artificial target in the 
upper left, can be resolved given. It is also likely that absolute changes in the optical 
properties can be better resolved than the absolute magnitude of the property itself, whereas 
only the latter case has been tested here. Further study is needed to determine the increase in 
image quality which is possible with incorporation of MRI-derived initial guesses of the 
optical property distribution, and other enhancments of the reconstruction process through 
adaptive spatial regularization or constraints [7, 40]. 

5, Conclusion 

In general, increasing the number of measurements improves the quality of the reconstructed 
images and the image accuracy. Care must be taken to ensure that the projections used 
provide a homogenous sampling of the entire tissue region of interest, and that the 
projections maximize the number of angles over which measurments are recorded. The 
ability to resolve objects improves as the distribution of angular projections increases through 
the sample; so, the maximum angular coverage of a tissue results from the circular geometry 
with liill sampling of the reflectance and transmittance light, which is demonstrated to have 
the best image quality in these simulations. The technological simplification of using a fan 
beam approach is attractive and comes at only a modest decrease in image quality, so this 
may provide the most practical method for optical imaging. The fan beam approach 
produces measurmements which are alt near the same order of magnitude in intensity so that 
variable attenuation is not required in the data acqusition system. 

A comparison of slab-geometry imaging to circular-geometry imaging with single 
objects does not demonstrate a major difference in the image quality, but the ability to image 
layered objects with the slab-geometry was seriously compromised. From these simulations. 
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we anticipate that the slab geometiy is not well suited to imaging layered structures, and this 
may be due to the relatively few projections at varying angles through the layers which exist. 
Sub-surface imaging produces the lowest quality images, but is able to quantify absorbing 
perturbations near the surface. Objects lower than 2 cm are not likely to be imaged 
accurately by subsurface methods unless the overlying surface is curved. Curved surfaces 
allow the light projections to probe deeper into the tissue and can provide accurate images of 
sections of a tissue. The image quality for all geometries presented in this paper is 
compromised by the generality of the inversion algorithm, so that some improvements can 
potentially be incorporated for each particular geometry through spatially distributed 
regularization methods or adaptive regularization [15], 
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In vivo local determination of tissue optical 
properties: applications to human brain 



Frederic Bevilacqua, Dominique Piguet, Pierre Marquet, Jeffrey D. Gross, 
Bruce J. Tromberg, and Christian Depeursinge 



Ivocal and superficial near-infrared (NIE) optical-property characterization of turbid biological tissues can 
be achieved by measurement of spatially resolved diffuse reflectance at small source-detector separa- 
tions (<1.4 mm). However, in these conditions the inverse problem, i.e., calculation of localized absorp- 
tion and the reduced scattering coefficients, is necessarily sensitive to the scattering phase function. 
This effect can be minimized if a new parameter of the phase function y, which depends on the first and 
the second moments of the phase function, is known. If 7 is unknown, an estimation of this parameter 
can be obtained by the measurement, but the uncertainty of the absorption coefficient is increased. A 
spatially resolved reflectance probe employing multiple detector fibers (0.3-1.4 mm from the source) is 
described. Monte Carlo simulations are used to determine y, the reduced scattering and absorption 
coefficients from reflectance data. Probe performance is assessed by measurements on phantoms, the 
optical properties of which were measured by other techniques [frequency domain photon migration 
(FDPM) and spatially resolved transmittancej. Our results show that changes in the absorption coef- 
ficient, the reduced scattering coefficient, and 7 can be measured to within ±0.005 mm '\ ±0.05 mm"^, 
and ±0.2, respectively. In vivo measurements performed intraoperatively on a human skull and brain 
are reported for four NIK wavelengths (674, 811, 849, 956 nm) when the spatially resolved probe and 
FDPM are used. The spatially resolved probe shows optimum measurement sensitivity in the mea- 
surement volume immediately beneath the probe (typically 1 mm^ in tissues), whereas. FDPM typically 
samples larger regions of tissues. Optical-property values for human skull, white matter, scar tissue, 
optic nerve, and tumors are reported that show distinct absorption and scattering differences between 
structures and a dependence on the pbase-fimction parameter 7, © 1999 Optical Society of America 
OCIS codes: 170.7050, 290.4210, 170.3660, 170.5280, 170.4550, 170.6510. 



1. Introduction 

Probing the optical properties of biological tissues has 
a major effect in several medical applications for di- 
agnosis and therapy. For example, the knowledge of 
these properties is necessary for optimizing tech- 
niques such as near-IR spectroscopy or photodynamic 
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therapy. The scattering and the absorption charac- 
teristics of many different kinds of tissues have been 
reported in the literature.^ However, they have 
been most often measured in vitro. Because of un- 
avoidable alterations in excised samples, such as 
blood drainage, structural alterations, and tempera- 
ture changes, these values are questionable and in 
vivo measurements are preferred. 

Measurement of optical properties performed in 
vivo could also be used as a diagnostic tool, which is 
complementary to other optical biopsy techniques, e.g., 
tissue autofluorescence. Indeed light scattering and 
absorption can provide information both on tissue 
structure and on chromophore content, features that 
can be used to distinguish between normal tissues, 
malignant lesions, and other pathologies. For exam- 
ple, hemoglobin and water content have been found to 
be significantly different in normal and cancerous tis- 
sues.2'3 Differentiation between normal and malig- 
nant bladder tissues* was fotmd to be pc^sible &om the 
elastic scattering and absorption properties. 
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Researchers have proposed different methods to 
determine quantitatively absorption and reduced 
scattering coefficients in vivo, by using spatially^'' 
and/or temporally resolved measurements. ^-^'^-^ 
Besides organs such as breast or neonatal brain, 
which can be transilluminated, measurements of 
thick tissues can be made in reflectance geometry. 
The case in which source- detector separations are 
larger than several transport mean-free paths, corre- 
sponding typically to distances greater than 5 mm for 
biologic^ tissues, has been extensively studied theo- 
retically and experimentally. Diffusion theory or 
Monte Carlo simulations have been commonly used 
to relate the measured light intensity to the optical 
coefficients. In particular, optimized source- 
detector separations have been calculated by differ- 
ent ,authors^'i° to ensure the best sensitivity to 
absorption and scattering coefficients from spatially 
resolved reflectance data. For typical turbid tissues , 
the optimal determination of both absoiption and 
scattering properties requires reflecteince measure- 
ment at small and large distances. For example, 
measurements were performed by Farrell et al.^ us- 
ing distances from 1 to 10 mm and by Bays et al.^ 
using distances from 3.6 to 15 mm or 2 to 14 mm. 

Nevertheless in all these studies the turbid me- 
dium is considered homogeneous or made of homoge- 
neous layers. Our approach is different. We wish 
to characterize optically a small volume of tissue, of 
the order of a few cubic millimeters, possibly distinct 
from surrounding tissues. Therefore our goal is to 
differentiate a small tissue heterogeneity instead of 
determining average optical properties of a large vol- 
ume of tissue. To achieve this, we chose to perform 
spatially resolved reflectance measurements with 
only small source- detector separations, from 0,3 to 
1.4 mm, even if these distances are not optimal for 
absorption determination. 

Mourant et a/.*'^^ have shown that the absorption 
coefficient can be estimated from measurement at a 
single distance of —1.7 mm, assuming the scattering 
coefficient to be in a certain range. They have also 
shown that measurements at a single shorter dis- 
tance (-0,2-0.3 mm) allows monitoring of spectro- 
scopic variations in the scattering properties of 
tissues. In this research we address a general case 
in which both tissue scattering and absorption prop- 
erties are simultaneously estunated from measure- 
ments at distances smaller than 1.4 mm. 

The theoretical framework for interpretation of the 
measured profile in terms of absorption and scatter- 
ing properties^2'^3 was given in a previous study 
based on Monte Carlo simxilations. The role of the 
phase function at short source- detector separations 
was carefully studied. In particular, we showed 
that, besides absorption and reduced scattering coef- 
ficients, a parameter depending on the first and the 
second moment of the phase function must be taken 
into account for source-detector separation ranging 
approximately from 0.3 to 10 transjKirt mean free 
paths. This analysis is more complete compared 
with previous research in which only the first mo- 



ment of the phase function (the anisotropy factor) 
was considered. 

First we describe the probe design and the Monte 
Carlo model developed to simulate the measured pro- 
file. Second we discuss the average tissue volume 
probed and sensitivity to medium boundaries from 
experiments and simulations. Thfrd we describe 
how the scattering and the absorption properties of 
tissues can be deduced from reflectance data. The 
accuracy of the proposed procedure is demonstrated 
on phantoms, the optical properties of which are mea- 
sured by others techniques. Finally we present and 
discuss measurements of human brain obtained in 
vivo. These clinical measurements were performed 
in parallel with a complementary method, frequency 
domain photon migration (FDPM),^'^ which probes a 
larger tissue volume. Optical properties from 
FDPM measurements are compared with values that 
we obtained by the spatially resolved method de- 
scribed here, using small source- detector separa- 
tions. 

2. Materials and Methods 

A. Definitions 

The spatially resolved reflectance is denoted J?(p), 
where p is the source- detector separation. It is de- 
fined by the backscattered power received by a detec- 
tor per unit area for a source of power xmity. In our 
measurements p ranges between 0.3 and 1.4 mm. 

The optical properties of tissues are the average 
refractive index n of the medium, the absorption co- 
efficient (JL„, the scattering coefficient [JL.„ and the 
phase function j)(cos 9), where 6 is the scattering 
angle, The phase function is the probability density 
function for cos 0. We consider the refractive index 
of tissuesi* as a constant, n = 1,4. 

It is also useful to define the reduced scattering 
coefficient fx^' - |ji,(l - g) and the transport mean 
free path, mfp' = (|x^' + it^o)~^, where^ is the anisot- 
ropy factor and is defined as the average of cos G. 
Generally the reduced scattering coefficient (x^' and 
the absorption coefficient are used to characterize 
optically thick tissue. Indeed for a high albedo me- 
dium the fight fluence rate depends only on jx^' and 
at distances of several transport mean fr^e paths 
(typically p > 5 mm for tissues) from the source (dif- 
fusion approximation). Therefore use of \x.J and 
is a natural choice if measurements are performed at 
such distances. Because we want to take measure- 
ments at closer distances, in the range of one trans- 
port mean free path we expect that some parameters 
of the phase function must be taken into account. 
This theoretical problem was fully studied with 
Monte Carlo simulations and reported in Refs. 12 and 
13. We consider below only the main implications of 
this research. 

For distances between 0,3 and 10 transport mean 
free paths we found that the reflectance curve de- 
pends on 11^, iXj,', and a third parameter 7 = (1 - 
^2)/(l ~ ^i)> where and g2 are, respectively, the 



4940 



APPLIED OPTICS / Vol. 38, No. 22 / 1 August 1999 



first and the second moments of the phase function. 
The n."^ moment g„ is defined as:i^'*^ 

gn = j PnCcos e)p(cos e)d(cos 9), (1) 

where P„ is the Legendre polynomial of order n. 
Note that the first moment is the anisotropy factor 
g. The parameter 7 is derived from the second-order 
similarity relations derived by Wyman et al.,^^ which 
are valid for a second-order anisotropic radiance (i.e., 
the radiance can be expanded as a series of Legendre 
poljmomials of the order of at most 2). For compar- 
ison, note that the diffusion approximation and first- 
order similarity relation correspond to a first-order 
anisotropic radiance (i.e., the radiance can be ex- 
panded as a series of Legendre polynomials of the 
order of at most 1). The role of the parameter, 7 = 
(1 — ^2)/(l gi): implies that the anisotropy factory 
(^gi) alone is not sufficient to predict correctly the 
reflectance curve close to the source. 

From this analysis, each tissue can be potentially 
characterized by three parameters: ix^,', and 7. 
However, owing to the restricted range of the source- 
detector distances that we want to use, the simulta- 
neous determination of \x„, ix,', and 7 is not always 
possible with a high degree of accuracy. The achiev- 
able accuracy depends on the optical properties them- 
selves and on experimental uncertainties. This 
problem is addressed in Section 3 for both phantom 
and tissue measurements. 

The parameter 7 may give interesting information 
about the tissue structure. Indeed, as measured by 
several authors,^''-^'' the tissue phase function can be 
approximated by summing a highly anisotropic phase 
fbnction, Pha(cos 6), due to large particles, plus a low 
anisotropic phase function, j?£a(cos 6), due to small 
particles: 

Pii»sue(cos e) = (1 - a)pHA(cos e) -1- api^(cos 9). (2) 

The coefficient ct is introduced to guarantee the nor- 
malization of Pt,isgue(COS e). 

The first term p^^(cos 6) has been fitted to the 
Rayleigh-Gans phase functional qj. iq ^j^g Henyey- 
Greenstein phase function,!'''.!^ j3j,jq(cos B): 

Pho(cos 9) = — \ — — ^ . (3) 

2(1+^hg'-2^hgC0s9P 

The moments of the Henyey-Green stein phase func- 
tion are given hy - ^hg" {n > Q). 

The second term, pi^f^{cos 9), can be interpreted as 
a Rayleigh-scattering contribution, valid for small 
scatterers compared with the wavelength. It has 
been approached by a purely isotropic term^'^'^^ or a 
Henyey-Greenstein phase function with a low nega- 
tive g value. 1^ We propose using the Rayleigh phase 
functioni^®: 

iJiuyioigh(cos9) = ^(l + coa^e). (4) 
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Fig. 1. Probe for measuring spatially resolved reflectance. 



The moments of the Rayleigh phase function aregi = 
0,g, = 0.1,^3 = 0, .... 

Rayleigh scattering contributes nothing to the first 
moment gi but only to the second moment gg- The 
moments of the tissue phase function ^^^^^^^(cos 9) 
given by Eq. (2) are therefore 

gi = (l-a)^HG, ^2 = (1 - + 0.1a. (5) 

Therefore we see that parameter 7 is influenced by 
the relative concentration of Rayleigh scatterer a, 
which should depend on tissue structure. 

Published phase-function data suggest possible 
values of 7 for biological tissues. I'^'-^o These phase 
functions have been measured in goniometric exper- 
iments on thin samples. Many artifacts can affect 
these measurements, such as tissue preparation and 
tissue thickness, and these results should be used 
with caution. The phase function reported by 
Jacques et al." for human dermis at 633 nm leads to 
7 = 1.4. The phase functions of white and gray mat- 
ter at \ = 750 nm measured by van der Zee et al.'^^ 
give values of —1.5. Therefore we performed Pmi-ig 
simulations with 7 = 1.5 as a starting point. As will 
be discussed, fitting our experimental data to simu- 
lations performed with different values of 7 estima- 
tion permits a reasonable parameter estimation, 

B. Experimental Setup 

The probe used for measuring the spatially resolved 
reflectance is described in Fig. 1. It is a linear array 
of optical fibers (core diameter, 200 |a,m; N.A. = 0.37 
in air). Two source fibers can be used to illuminate 
the tissue. They are disposed symmetrically vrith 
regard to the collecting fibers. If the sample is ho- 
mogeneous, the reflectance curve is identical with 
either illuminating fiber. Therefore comparing the 
two curves tests the heterogeneity of the investigated 
tissue region or detects obstructions beneath the il- 
luminating fibers. If the two curves are close (typi- 
cally with differences of less than 10%), the 
measurement is validated and the average of the two 
curves is calculated. 

The illuminating fibers are slid into small 
stainless-steel tubes to avoid direct light coupling 
with the collecting fibers. The coupling between 
each collecting fiber has been experimentally mea- 
siu'ed and foxmd to be less than 2%. The fiber array 
is set in a stainless-steel tube 2.5 mm in diameter and 
20 cm long. The tube is filled with an optically clear 
adhesive (black adhesive may also be used). The 
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Fig. 2, Experimental setup. 



probe is rigid, which allows for easier handling by the 
physician, during surgery, for example. The whole 
probe can be sterilized. The experimental setup is 
shown in Fig. 2, An optical switch (Dicon Model 
GP700) is used to select the illuminating fiber from 
different sources . For the brain measurements, four 
laser diodes emitting at 674, 811, 849, and 956 rnn 
were used (SDL, Inc. Models 7421, 5420, 5421, and 
6321, respectively). Two other laser diodes were 
used for phantom measurements, emitting at 675 
and 828 nm (ILEE LDA 2011 and 1805, respectively). 
The six fibers that collect the backscattered light are 
imaged on a linear CCD (Hamamatsu S3921). The 
signal is digitized by a 12-bit analog-to-digital card. 
Only one measurement, which takes —0.1 s, is then 
needed to measure simultaneously the intensity col- 
lected by the six fibers. The entire system is con- 
trolled by a personal computer. 

Transmission differences between each fiber are 
corrected by using a measurement on a turbid phan- 
tom illuminated uniformly. Immediately after each 
reflectance measurement, the background light is au- 
tomatically measured and then subtracted from the 
reflectance signal. To minimize the background 
light, a long-pass filter (X > 650 nm) is placed be- 
tween the end of the bundle and the CCD. Even 
during open surgery where ambient light is substan- 
tial, the measured background ranged between 1% 
and 10% at the farthest fibers, depending on the 
wavelength used. 

C. Monte Carlo Simulations 

A model of photon migration in tissues is necessary to 
define the relationship between the measured reflec- 
tance and the optical properties. Analytical solu- 
tions firom the diffusion equation are not appropriate 
in our case because we are interested in the reflec- 
tance close to the source, at a distance comparable 
with the transport mean ii-ee path (mlp').^ '' We per- 
formed Monte Carlo simulations to predict the mea- 
sured reflectance of an homogeneous semi-infinite 
turbid media. The code that we used was exten- 
sively tested. ^^'^^-^^ Any phase fiinction can be im- 
plemented in discretized form. 

Our simulations take into account the exact diam- 
eter of the illuminating and the collecting fibers as 
well as their numerical apertures (N.A. = 0,26 in 
tissue). However, the distortion of the signal due to 
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Fig. 3, Comparison between the exact and the simplified bound- 
ary conditions. The optical coefficients are fi/ = 1 mm"\ = 
0.01 mni"\ and-y = 1.9. 



the size of the fibers that we used (0 = 200 nm) has 
an almost negligible influence on the reflectance curve. 
The mismatch of the index of refraction at the sur- 
face of the medium is also taken into account in our 
simulation by using the Fresnel law for each photon 
reaching the surface. Some simxilations have been 
performed with the exact geometry of the probe as 
described in Fig, 1, taking into account the mismatch 
of the index of refi-action between the probe adhesive 
{n = 1.5) and the sample as well as the mismatch of 
the index of refraction between the air and the sam- 
ple outside the probe (p > 1.8 mm). Photons re- 
fracted inside the limien of the stainless-steel tube 
are considered absorbed. In Fig. 3 the reflectance 
obtained with the exact probe configuration is com- 
pared with the simplified case of the semi-infinite 



e(«m. 



= 1.5). The exact-probe- 



boundary condition is extrapolated beyond the sur- 
face of the probe to where the boundary is tissue/ air. 
As expected, there are important differences between 
the exact and the simplified cases for p > 1.5 mm, 
close to the limit of the probe. The decreased reflec- 
tance at p > 1.8 for the exact case is due to the 
increase in internal reflection at the interface be- 
tween the medium and the air. However, the differ- 
ences between the exact and the simplified cases are 
less than a few percent for p < 1.4 mm, corresponding 
to the region where the measurements are per- 
formed. 

The simplified case, assuming cylindrical symme- 
try, is computationally much less tune-consuming 
than the exact configuration. Exact-configuration 
.simulations require -- 10 times the number of photons 
for statistical errors to be achieved comparable for 
the simplified case. Indeed in the simplified case the 
cylindrical symmetry allows for the use of annular 
detectors to compute the reflectance. The advan- 
tage of annular detectors is their larger area, com- 
pared with the detector area used in the 
configuration of the exact-boundary condition config- 
uration. Following the result illustrated in Fig, 3, 
we decided to employ only the simplified semi-infinite 
condition and be restricted to p < 1,4 mm. 
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Fig. 4, Comparison between measurements in semi-inflnite and 
in inimite media. The medium is Intralipid. The optical coeffi- 
cients were measured by the FDPM technique: n,' = 1.2 mm"', 
(JL„ = 0.0005 mm-' (X = 675 nm). 



3. Results and Discussion 

A. Boundary Effects 

During in vivo investigations the ideal case of a me- 
dium with a perfect plane surface is never realized. 
Moreover the probe could be slightly pushed inside 
the tissue. Therefore the effect of the medium 
boundary on the measured reflectance is important to 
quantify experimentally. In Fig. 4 we show two 
measurements in Intralipid (\ = 675 nm), one at the 
surface and the second inside the medium. The dif- 
ference between these two measurements is less than 
5%, which could seem very surprising at first glance. 
However, as shown by simulation in Fig. 3, the 
boundary condition outside the probe has only a weak 
effect on the intensity measured by the six fibers. 
The boundary condition created by the probe itself, 
i.e., the index of refraction mismatch between the 
medium and the adhesive inside the probe, is much 
more critical. 

The negligible effect of the tissue boundaries 
shown in this section is an important advantage for 
clinical investigation. This experiment also clearly 
demonstrates that the sample volume investigated is 
principally confined to the region just beneath the 
probe surface. This point is developed further in 
Subsection 3.B. 

B. Depth of Tissue Investigated 

To quantify more precisely the volume of tissue 
probed in our technique, we determined the depth 
below the surface of each scattering event in the sim- 
ulation. With this information we determined the 
average depth of all the scattering events for each 
detected photon, which we present as a probability 
density function in Fig. 5, To permit more general 
statements, distances are expressed here in mfp' 
units. For typical turbid tissue^ and near-IR wave- 
lengths, |jLs' is around 1 mm"-^ and is less than 0.2 
mm~^, which means that 1 mfp' =« 1 mm. Figure 5 
shows that the average depth of scattering is around 
0.8 mfjp'. Moreover it shows that for typical tissue 



per delected photon [mfp'] 

Fig. 5, Probability density of the mean depth of scatterii^ events 

per detected photon for the case of 7 = 1.9. 



optical properties, structures located beneath the 
probe at a depth greater than 2 mfp' are not likely to 
contribute significantly to the measured signal (for 
p < 1.5 mfp'). These results are consistent with the 
average depth computed by Weiss et al.^ using 
Monte Carlo simulations. However, note that their 
simulations do not take into account any refractive- 
index mismatch or a restricted nimierical aperture of 
the detector. 

To evaluate the maximum depth of photons paths, 
we performed experiments in Intralipid, placing the 
probe at the liquid surface and moving an absorbing 
plate placed horizontally, as described in Fig, 6. The 
reflectance J?(p, d) was measured for varjring thick- 
ness d. The ratio iJ(p, d)/R{p, d = «) is reported in 
Fig. 6. Figure 6 shows that the intensity of the re- 
flectance is decreased by —20% if the medium is 2 
mfp' thick and by ~10% if it is 3 mfp' thick {for the 
albedo and the range of p considered) . These results 
imply that 80% or 90% of the photons do not pene- 
trate deeper than 2 or 3 mfp', respectively, into the 
medium. This experiment complements our simu- 
lation result (Fig. 5) where the aver^ depth of scat- 




Fig. 6. Effect of the Intralipid thickness on the reflectance. The 
setup for the investigation of Intralipid with varying thickness is 
illustrated on the right-hand side. The optical coefficients of In- 
tralipid were measured by the FDPM technique (\ = 956 nm). 
The reduced albedo is a' = m,,7((j.,' + |xj = 0.98. 
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Fig. 7. Calibration measurement of a siloxane phantom and mea- 
surement test of a microsphere suspension. The siloxane phan- 
tom optical piroperties were measured with PDPM; = 
0.00024 ± 0.00002 mm-^ ft.' = 1-82 ± 0.01 mm"^ (X = 674 nm). 
The siloxane phantom measurement is multiplied by the calibra- 
tion fa<:tor to fit the corresponding simxJation (7 = 1.8). The 
measurement on microsphere suspension is multiplied by the cal- 
ibration factor derived from the siloxane phantom measurement. 
The simulation, corresponding to the microsphere suspension mea- 
surement, is performed by using the optical properties derived 
from Mie theory; (jl^' = 1.0 mm~^, = 0.0005 mm~^, Mie phase 
function (gi = 0.916, y = 2.2). 



tering was estimated to be ~1 mfp'. Thus for typical 
biological tissues our measurements are maioly sen- 
sitive to the region of tissue located within 2 mm of 
the surface and the investigated volume is of the 
order of 1 mm^. 

C. Calibration and Test on IVIicrosphere Suspension 
To perform absolute intensity measurements, cali- 
bration is performed on a solid turbid siloxane phan- 
tom of known fx^ and jjl^' (determined independently 
by FDPM).2.3 Value -y = 1.8 ± 0.1 was determined 
for the calibration phantom by consideration of the 
value that gives the best fit to the experimental 
curve, |a,„ and [x.,' being fixed (\ — 674-956 nm). 

As shown in Fig. 7 the phantom measurement, 
scaled by a single factor applied equally to all fibers, 
fits well the Monte Carlo simulation performed with 
phantom coefficients. This factor is defined as the 
calibration factor for a given wavelength. The cali- 
bration was performed at the end of each set of clin- 
ical measurements. 

Experiments on microsphere suspensions (polysty- 
rene sphere 0 1.072 ± 0,019 |im) were performed to 
assess the accuracy of our theoretical model and the 
calibration method. The scattering coefficient and 
the phase function of such turbid media can be known 
precisely by using Mie theory.^s Because no dye was 
added to the suspension, was considered to be 
equal to the water absorption. In Fig. 7 measure- 
ment of the reflectance is compared with a simulation 
computed with the microsphere suspension coeffi- 
cients ((jl/ = 1.0 mm""\ |JL„ = 0.00041 mm"\ 7 = 2.2), 
The excellent agreement found here between exper- 
iments and simulation confirms the accuracy of our 
simulations and the validity of our calibration proce- 
dure. 
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Fig. 8. Relationship between parameters ij(p = 1 mm) and la^ In 
i? (p = 1 mm)| and the optica) coefficients fi,' and m,„ for cases <s£y ~ 
1.5 and 1.9. 



D. inverse Problem 

Our goal is to solve the inverse problem that consists 
of extracting optical coefficients from the reflectance 
data. Measurements of the reflectance intensity 
R{p) and the slope of In R{p) [denoted d In Rip)], 
determined at a distance p = 1 mm, can be used to 
derive p,/ and for a given 7 value. Figure 8 shows 
graphically the relationship between ^j.^', and p^ and 
the two parameters, R{p = 1 mm) and ja,, In J?(p = 1 
mm)|. To illustrate the influence of parameter 7, 
two examples are superimposed: 7 = 1.5 and 7 = 
1.9. We see clearly in Fig. 8 that p^', and p,„ cannot 
be determined uniquely from the two parameters, 
R{p = 1 mm) and |9p In R{p = 1 mm)|, if 7 is unknown. 
This indetermination may be resolved by the values 
of R{p) and/or In R{p)\ at other distances. There- 
fore the following procedure was defined for tissue 
measurements: 

(1) Determination of p.,' and Po from J?(p = 1 mm) 
and In i?(p = 1 mm)| for discrete values of 7 (for 
example, 7 = 1,0, 1.5, 1.9, 2.2), 

(2) Simulations with different sets of p^' and p„ 
obtained from Eq. (1), 

(3) Comparison between the simulations and the 
reflectance profile for distances 0.35 < p < 1.4 mm. 

This last step permits us to determine the value of 
7 that gives the best fit. Points (l)-(3) can be done 
iteratively to evaluate 7 more precisely. The preci- 
sion that can be obtained depends on the optical co- 
efficients themselves and on the experimental 
uncertainties. Nevertheless two important techni- 
cal points should be noted: 

First, determination of p^,' is only weakly influ- 
enced by 7 for p^' close to 1 mm"^. Indeed in Fig. 8 
the differences induced by 7 = 1.5 or 7 = 1.9 are 
typically 10% for (jl^'. In contrast, the absolute de- 
termination of p„ is critically sensitive to 7. How- 
ever, if 7 remains constant, relative variations of p„ 
can still be precisely evaluated. This point is dis- 
cussed with the results obtained from Intralipid 
measurements in Subsection 3.E, Metabolism moni- 
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toring or drug monitoring could therefore be a poten- 
tial application of such a probe. 

Second, the experimental determination of In 
Rip = 1 mm)| requires measurements at different 
distances close to 1 mm. To minimize errors on In 
i?(p = 1 mm)| due to experimental artifacts, we per- 
formed a fit of the reflectance curi-e (0.5 mm < p 1.4 
mm) with the function m-ip'"' expf-mgp), which was 
always found to fit Monte Carlo simulations well for 
this restricted range of distances. (The same func- 
tion was also proposed by Bolt and ten Bosch.*^) Pa- 
rameters R{p = 1 mm) and \dp In R{p = 1 mm)| are 
then derived from the fit. Once optical coefficients 
|x„ and are derived, the validity of this procedure 
is double-checked by comparing the curve obtained 
from the Monte Carlo simulation with the experimen- 
tal profile. 

E. Phantom Measurements 

We present in this subsection measurements on tis- 
suelike phantoms (Intralipid and imcrosphere sus- 
pensions). Measurements of phantoms with 
varying jjl^ and |Xj.' values were used to test the in- 
version procedure. For the Intralipid measure- 
ments we considered the y value as a priori not 
known, even though 7 can be estimated from pub- 
lished research (described below). The phantoms 
made of Intralipid and dye [iron(III) ferrocyanide] at 
different concentrations were calibrated by the 
FDPM technique (performed at large source- 
detector separations and therefore insensitive to 7), 
For Intralipid measurements, reflectance profiles 
were found to fit best when the parameter 7 was 
between 1.6 and 1.8 for \ = 674-849 nm. The 
values derived from our method are plotted in Fig. 
9(a) versus the values obtained by the FDPM tech- 
nique. Figure 9(a) shows that 7 = 1.8 leads to an 
overestimation of by -0.005 mm"-', whereas 7 = 

1.6 leads to an underestimation of tJL„ by —0.01 mm~^. 
This example clearly illustrates that determination 
of the absolute values of fL„ are sensitive to 7. How- 
ever, changes in the absorption coefficient can be 
measured to within ±0.005 mm"'. 

The accuracy of the inversion procedure on pi^.' is 
illustrated in Fig. 9(b). In the case of jx.,', as men- 
tioned above, the influence of 7 is weaker than in the 
case of |JL„. A difference of only 5% is found on the p,^ ' 
if 7 = 1.6 is used instead of 1.8. This is also approx- 
imately the variation in |jLj' values obtained by mul- 
tiple FDPM measurements. We also found that in 
the case of a constant Intralipid concentration the 
measured values of vary by less than 2% when |ui,„ 
is increased from 0.012 to 0.5 mm~^ by the addition of 
dye, which proves that the data-inversion procedure 
effectively uncouples p.^ and p,^'. 

Prom our measurements, 7 can be estimated to 7 = 

1.7 ± 0.1 (\ = 675 nm). The calculation when Mie 
theory was used and the size distribution given by 
van Staveren et alJ' leads to 7 = 1.89 (\ = 675 nm). 
These values are in good agreement considering that 
the actual size distribution of the sample may slightly 
vary from the sample measured by van Staveren et 
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Kg. 9. Comparison between (a) |x„, and (b) |x,' obtained hfy the 
FDPM technique and by probe measurements. Measurements 
are on Intralipid [{a) with dye; (b) without dye]. 



al. Note also that they found a similar discrepancy 
for between Mie theory and experiments. 

Further assessments of the inversion pTOcedure 
were performed on monodisperse microsphere sus- 
pensions. In this case we considered the 7 value as 
a priori known from Mie theory; 7 = 2.2 (polysty- 
rene, 0 = 1.053 ± 0.010 |xm, \ = 675 nm). The 
optical coefficients derived from our local reflectance 
measiu-ements are compared with values obtained 
with a spatially resolved transmittance method de- 
scribed elsewhere,^^'^'^ 

Figures 10(a) and 10(b) show and n,,' values 
obtained by both methods: spatially resolved trans- 
mittance and reflectance. Phantoms made of micro- 
sphere suspensions and ink at different 
concentrations were used. An excellent correlation 
is found between the p.„ and \lJ values obtained by 
the two methods. Small systematic differences (typ- 
ically 10%) are found when comparing absolute val- 
ues. They are mainly due to errors occurring in the 
calibration procedures (for both methods). A further 
comparison with Mie theory also showed a typical 
deviation of 2% for (jl^' (a maximum deviation of 10%). 
The relatively small errors that we foimd here are 
typical of errors fovmd when different techniques for 
measuring turbid-media optical properties are com- 
pared.™ Such errors could be avoided by multiple 
calibrations on several turbid samples of different 
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Fig. 10. Comparison between (a) and (b) p.,', obtained from 
spatially resolved transmittance (the method described in Refs. 22 
and 23) and from probe measurements. Measurements are on 
microsphere suspensions [(a) with ink; (b) without ink]. 



known optical properties. However, the accuracy of 
tissue measurements is subject to other major limi- 
tations due to, for example, their structure and het- 
erogeneity.29 Therefore more accurate calibration is 
not necessarily required for absolute tissue measure- 
ments, since we are mainly interested in observing 
optical-property diflferences between tissue physio- 
logical states, 

F. In vivo Measurements on Brain Tissues 
Clinical measurements of normal and malignant 
neiural tissues were recorded in vivo during brain 
surgery, ^'^ Two different cases are reported here. 
Case 1 was a 3-year-old male and case 2 was an 
8-year-old male. Different types of tissues were in- 
vestigated in each case. Several measurements 
(typically six) were always performed successively at 
a given location. The intensity fluctuations (typical- 
ly of the order of 10%) for these measurements were 
mainly due to tissue heterogeneity and slight probe 
movements. The average reflectance was calculated 
for each location as well as the standard deviation, 
Note that the uncertainty due to the apparatus, es- 
timated from measurements on a phantom, is much 



lower (<5%). Before each set d' measurements the 
blood from the surgical site was carefully irrigated 
away with saline and the probe was cleaned with a 
saline-damped sponge. The measurements pre- 
sented here were performed in parallel with 
frequency-domain measurements^'^ (FDPM) by using 
a source- detector separation of 10-14 mm. 

In case 1 measurements were performed on normeil 
cerebral cortex (frontal lobe and temporal lobe), optic 
nerve astrocytoma (size, =«1.3 cm), and normal optic 
nerve. In case 2 measurements were performed on 
the skull, deep cerebellar white matter with scar tis- 
sue (from a previous surgery), meduUoblastoma (size, 
=«3.8 cm), and deep cerebellar white matter. Tumor 
dimensions were estimated from conventional imag- 
ing techniques (i.e., computed tomography and/or 
magnetic resonance imaging). 

As discussed in Subsection 3.B, the depth probed is 
less than ~2 mm. For each type of tissue that we 
investigated the influence of surrounding tissues on 
the measurement is weak. In particular, only gray 
matter is investigated during measurements of the 
cerebral cortical surface. 

Figures 11(a) and 11(b) show the measured param- 
eters R{p = 1 mm) as a function of \d In R(p = 1 mm)| 
obtained for cases 1 and 2, respectively. They are 
similar to Fig. 8 except that the relationship between 
parameters i?(p = 1 mm) and \dp \nR{p = 1 mm)| and 
optical coefficients and ii^' is indicated only qual- 
itatively by two arrows for better clarity. (Different 
grids corresponding to different y values should be 
superimposed.) Quantitative results are reported in 
Tables 1 and 2. 

Figures 11(a) and 11(b) show that parameters 
Rip - 1 mm) and |3p \nR{p = 1 mm)| provide excellent 
discrimination between tissue types. The spectro- 
scopic signature on R{p = 1 mm) and |3p In R{p = 1 
mm) I should also be noted. The cortex and the skull 
exhibit less significant spectroscopic differences com- 
pared with tumor tissues such as the astrocytoma 
and the medulloblastoma. Thus parameters J?(p = 
1 mm) and la,, In E(p = 1 mm)| could be useful for 
optical biopsy. Nevertheless we believe that to ex- 
ploit these results fully the differences found must be 
explained in terms of scattering and absorption pa- 
rameters. These factors in turn can be used to un- 
derstand physiological and structural variations. 
The procedure described in Subsection 3.D was used 
to determine coefficients ix^', i-l^j, and 7 from the mea- 
sured curves. This procedure is fully illustrated for 
measurements of the temporal lobe. Results of op- 
tical coefficient calculations are summarized for all 
tissues in Tables 1 and 2. Values obtained in par- 
allel by the FDPM technique are also indicated. 

Figure 12 shows that the best fit to the cortex 
(temporal lobe) data is obtained with 7 = 1.9. 
Lower values, such as 7 = 1.5, led to impossible 
values for jR(p = 1 mm) and lnjR(p ^ 1 mm)| and 
should therefore be rejected. Larger values, such 
as 7 = 2.2, do not fit the reflectance data for dis- 
tances p < 0.8 mm. The 7 value is then estimated 
to be 7 = 1.9 ± 0.2. Note that 7 = 2.2 would lead 
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Fig. 11. Clinical measurements in vivo on human brain for (a) case 
1 and (b) case 2. Plot of fl(p = 1) and In ft(p = 1 mm)| for different 
types of brain tissue: (a) normal cortex (frontal and temporal lobe), 
optic nerve astrocytoma, and normal optic nerve; Go) skull, deep 
cerebellar white matter with scar tissue, medulloblastoma, and deep 
cerebellar white matter (surrounding normal tissue). 



to almost identical (jl^' (differences of less than 5%) 
and an overestimation of fjL„ of —0.02 mm"-'-. Tak- 
ing into account the uncertainties of the measure- 
ments [-5% on R{p = 1 mm) and In R{p = 1 mm)| 
for this tissue] and the uncertainty of 7 (±0.2), the 



error of the absolute value of iJi,^' is estimated to be 
5%. 

In some cases, as shown, for example, in Fig. 12, 
the nearest fiber does not perfectly fit the simulation. 
Different reasons may be suggested. First, this fiber 
is at the limit of the region where only the parameter 
7 must be taken into account. Therefore higher mo- 
ments of the phase function may also be considered. 
Second, this fiber probes the tissue more superficially 
than the other fibers (see Fig, 5), and the differences 
that occur at the nearest fiber may be due also to 
tissue heterogeneity. 

As reported in Table 1, with the exception of the 
normal optic nerve, the absorption coefficients at 811 
and 849 nm are lower than those obtained at\ = 674 
nm or \ = 956 nm. This result is consistent -with the 
fact that the main near-IR tissue chromophores, he- 
moglobin and water, have an absorption maximum at 
approximately X < 700 nm and X = 970 nm, respec- 
tively. ''i The overall absorption is much higher in 
the tumor than the cortex presumably because of the 
greater hemoglobin in the tumor. Tumors generally 
grow more new blood vessels and thus have a higher 
hemoglobin concentration. 

The jjL,^' of the tumor is similar to the cortex at \ = 
956 nm, However, the variation of jAj,' between X = 
674 nm and k = 956 nm is much larger for the tumor 
than for the cortex: AiJ.^' « 0,57 mm~^ for tumor 
and A[i^' 0, 10 mm^^ for the cortex. Such spectro- 
scopic variations may be attributable to structural 
differences between tissue types. Indeed such dif- 
ferences may depend on the average size or size dis- 
tribution of scattering structures within or between 
cells. 

Higher scattering is found for the normal optic 
nerve compared with all the other tissues in case 1. 
This is likely due to the presence of myelin. In vitro 
measurements have also shown that myelin contain- 
ing white matter exhibits a higher scattering coeffi- 



Tabte 1. Optical Properties of Norma! and Malignant Hi 



in Brain Tissue, Case 1 



Probe Measurements 



Type of Tissue 


'ivavejengin 
(nm) 




!^;( 






(mm-^) 


(mm ^) 


Cortex (frontal lobe) 


674 


1.9 ± 0.2 


1.00 


± 0.05 


<0.02 ± 0,01 


1.12 


0.0173 




811 


1.9 ± 0.2 


0.91 


± 0.05 


<0.01 ± 0.01 


0.74 


0.0182 




849 


1.9 ± 0.2 


0.92 


±0.05 


<0.01 ± 0.01 


0.74 


0.0185 




956 


1,9 ± 0.2 


0.89 


±0.05 


0,015 ± 0,01 


0.80 


0.0206 


Cortex (temporal lobe) 


674 


1,9 ± 0.2 


1.00 


±0.05 


0,02 ± 0,01 


0.99 


0,0179 




811 


1,9 ± 0.2 


0.82 


±0.05 


0,02 ± 0,01 


0.48 


0,0190 




849 


1,9 ± 0.2 




±0.05 


<0,01 ± 0,01 


0.45 


0,0179 




956 


1.9 ± 0.2 


0.82 


±0.05 


0.025 ± 0.01 


0.42 


0.0218 


Astrocytoma of optic nerve 


674 


1.7 ± 0.2 


1.25 


±0.10 


0.14 ± 0.03 


0.92 


0.0165 




811 


1.7 ± 0.2 


0.95 


± 0.10 


0.12 ± 0.03 


0.55 


0.0190 




849 


1.7 ± 0.2 


0.76 


±0.10 


0.09 ± 0.03 


0.59 


0.0191 




956 


1.7 ± 0.2 


0.73 


±0.10 


0.15 ± 0.03 


0.58 


0.0323 


Normal optic nerve 


674 


1.7 ± 0.2 


1.75 


± 0.20 


0.06 ± 0.03 


N/A 


N/A 




811 


1.7 ± 0.2 


N/A 


N/A 


N/A 


N/A 




849 


1.7 ± 0.2 


1.60 


±0.20 


0.08 ± 0.03 


N/A 


N/A 




956 


1,7 ± 0,2 


1,52 


±0,20 


0.07 ± 0.03 


N/A 


N/A 



'The unceartainty of the FDPM values, given by the fitting procedure (see Eefe. 2 and 3), is typically 2%. 
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Table 2. OpUcal Properties of Normal and Mafignant Human Brain Tissue, Case 2 



Probe Measurements FDPM* 

Wavelength — 

Type of Tissues (nro) 7 |a/ (mm"') (a^ (mm"') |x,' (mm"") (mm"') 



Skull 


674 


1.9 ± 0.2 


0.9 ± 0.1 


0.05 + 0.02 


1.19 


0.0208 




849 


1.9 ± 0.2 


0.9 ± 0,1 


0.05 ± 0.02 


0.91 


0.0215 




956 


1.9 ± 0.2 


0.85 ± 0,1 


0.05 ± 0.02 


0.77 


0.035S 


Cerebellar wbite matter 


674 


1.9 ± 0.2 


1.35 + 0.1 


0.25 ± 0.05 


1,34 


0.0165 




S49 


1.9 ± 0.2 


0.S5 ±0.1 


0.095 ± 0.02 


0.98 


0,0132 




956 


1.9 ± 0.2 


0.78 ± 0.1 


0.090 ± 0.02 


0.84 


0.0299 


Medtdloblastoma 


674 


1.9 ± 0.2 


1.40 ± 0.1 


0.26 ± 0.05 


1.05 


0.0120 




849 


1.9 ± 0.2 


1.07 ± 0.1 


0.10 ± 0,02 


0,66 


0.0079 




956 


1.9 ± 0.2 


0,4 ± 0.1 


0.075 ± 0.02 


0.54 


0.0239 


Cerebellar white matter with scar tissv 


les 674 


2.2 ± 0.2 


0.65 ± 0,05 


<0.02 


N/A 


N/A 




849 


2.2 ± 0.2 


0,80 ± 0,05 


<0.02 


N/A 


N/A 




956 


2.2 ± 0.2 


0.65 ± 0.05 


<0.02 


N/A 


N/A 



"The uncertainty of the FDPM values, given by the fitting procedure (see Hefe. 2 and 3), is typically 2.5%. 



cient compared with other tissue. ^^'i^'^^ The 
accuracy of the measurement on the optic nerve tis- 
sue can be affected by the high anisotropy of this type 
of tissue. Indeed it has been reported that hght 
propagation depends on whether the direction con- 
sidered is parallel or perpendicular to the nerve 
fibers.32-34 Interestingly, comparison between mea- 
surements with the two symmetric sources reveals 
larger differences than for the other tissues. This 
variation is represented in Fig. 11(a) by the relatively 
large uncertainties associated with the optic nerve 
measurements. 

In case 2 the first measurements were acquired 
directly from the skull. The difference in reduced 
scattering properties between X - 674 and 956 nm is 
very small (ix^.' = 0.9 and 0.85 mm~^, respectively) as 
shown in Fig. 11(b). Previously determined in vitro 
values^^ for pig skull were approximately twice as 
large. Such differences may be due to variations in 
water content and sample preparation between in 
vivo and in vitro studies. 

Case 2 data were also obtained from in series deep 
cerebellar white matter with scar tissue, tumor (me- 
diilloblastoma), and cerebellar white matter in the 
excavated tumor bed. Figure 11(b) shows that scar 




Fig. 12. Comparison between the spatially resolved reflectance 
curve measured in vivo on normal cortex (temporal lobe) and sim- 
ulations. 



tissue is well differentiated from other structures and 
is characterized by a high 7 value (7 = 2.2) and low 
absorption ((i,„ < 0.02 mm~-^), which is consistent 
with the low vascularization of such tissue. In con- 
trast, larger la,^, smaller 7, and large spectroscopic 
differences in are foimd in white matter and me- 
dulloblastoma. As in case 1 the large [x^ values ob- 
tained for the tumor can be related to higher 
hemoglobin content, as is often found for cancerous 
tissues. Taking into account measurement variabil- 
ity, no signiiicant differences are found between the 
tumor and the surrounding normal white matter. 
However, it is not clear that the so-called normal 
tissue measured at the surgical boarder could be con- 
sidered disease-free and unaffected by the tumor vas- 
culature. 

Some of the different types of tissue examined here 
showed a d(;pendonce on the phase-function param- 
eter 7. This confirms that 7 may be a valuable pa- 
rameter for tissue characterization. Nevertheless 
an improvement in the accuracy of the 7 determina- 
tion would be necessary to conclude the importance of 
this parameter. 

The optical coefficients that we obtained can be 
compareil with measurements performed simulta- 
neously with the FDPM technique. One should 
keep in mind that the depth investigated by the 
FDPM technique is greater compared with the spa- 
tially resolved technique described here. Generally 
the p,^.' values are similar between the two methods, 
whereas more differences are found for fx^. For the 
cortex and the skull measurements both ix,J and 
values are generally in good agreement. The de- 
crease in 10,^' from \ = 674-956 nm is more pro- 
nounced in the FDPM data. In contrast, important 
differences are found for results where tumor val- 
ues (astrocytoma and meduUoblastoma) obtained by 
FDPM are significantly lower. This can be ex- 
plained by the sensitivity of the spatially resolved 
probe to liie high local hemoglobin content that can 
be resolved only by the small source- detector sepa- 
rations. La contrast the large source-detector sep- 
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aration employed by the FDPM probe interrogates 
much greater larger-tissue volumes and hence mea- 
sures average optical properties firom multiple struc- 
tures (e.g., normal + malignant). 

4. Conclusions 

The purpose of this research was to assess the per- 
formance of spatially resolved reflectance by using 
short source- detector separations (<1.4 nmi). 
Monte Carlo simulations, the accuracy of which were 
confirmed by experiments on tissue phantoms, were 
used to establish the correspondence between the 
measured reflectance and the optical properties. 

Optical properties that can be determined by this 
technique are the absorption coefficient p.^' 
duced scattering coefficient ji^', and a parameter of 
the phase function 7 = (1 -g^/\l -^i), where 
and ^2 *re the first and the second moment of the 
phase function. Experiments on calibrated In- 
tralipid solutions and microsphere suspensions 
showed that changes in the absorption coefficient, 
reduced scattering coefdcient, and 7 can be measured 
to within ±0.005 mm~^, ±0.05 mm"^, and ±0.2 re- 
spectively. Systematic errors are possible if the pa- 
rameter 7 is not determined with sufficient accuracy. 
These performances could be improved by lowering 
the uncertainty on reflectance measurements. 

Experiments and simulations helped to define the 
average volume probed by this technique. For typ- 
ical tissues the average probe depth is typically --1 
mm and the influence of layers located below 3 mm is 
negligible. 

Finally in vivo measurements on human brain 
showed that excellent discrimination can be obtained 
between different types of neural tissues, normal and 
abnormal. Gfood correlation has been found be- 
tween spatially resolved reflectance and simulta- 
neous measurements performed by FDPM. These 
two techniques offer interesting complementary fea- 
tures. The spatially resolved probe can potentially 
provide better differentiation between different types 
of tissue owing to its sensitivity to local structure, 
because substantially a smaller volume of tissue is 
probed. On the other hand, owing to physical limi- 
tations imposed by large near-IR mean absorption 
lengths in tissue, the precision for the ix^ estimate is 
likely to be worse. Consequently, the short distance, 
spatially resolved technique appears to be well suited 
for clinical settings that require rapid localized tissue 
identification, such as endoscopic or a needle-based 
optical biopsy, and intraoperative tissue mapping for 
surgical guidance. 
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Near-infrared (NIR) spectroscopy has been used in 
studies of the cerebral hemodynamic response to vi- 
sual processing. In this paper, we present theoretical 
results from finite element and Monte Carlo modeling 
in order to help understand the contribution to the 
NIR signal from different parts of the head. The results 
from the models show that at the typical optode spac- 
ings used in these studies, an infrared spectroscopy 
measurement of intensity is sensitive to the outer 1-2 
mm of the cortical gray matter and the partial optical 
path lei^th in the gray matter is approximately 10 mm, 
compared with a total optical path length of 400 mm. 
When the NIR measurement is of change in mean 
photon arrival time (or phase shift), the signal comes 
from the upper 2-4 mm of the cortical surface and 
there is an increased lateral spread of the contributing 
tissue. We predict that for a 4-cm separation of input 
and detection optodes at 800 nm, a 1 |i.M change in 
hemoglobin concentration in the cortex corresponds 
to an attenuation change of approximately 0.001 OD 
(optical density) or 1 ps mean time change. Movement 
of the brain caused by this increase in volume will 
cause an absorption change of approximately half this 
magnitude, but does not affect the photon arrival time 
at 4-cm spacing. A discrepancy between the predicted 
and the experimentally measured intensities may sup- 
port the supposition that the NIR signed is actually 
very sensitive to changes occurring in the pial cerebral 
vessels lying on the brain surface. isss Academic Press 



INTRODUCTION 

Near-infrared (NIR) spectroscopy can be used to 
noninvasively monitor variations in cerebral blood oxy- 
genation by measuring changes in the attenuation of 
NIR light passing through tissue (Jobsis, 1977; Ed- 

' To whom correspondence should be addressed. 



wards et al., 1988; Hampson and Piantadosi, 1988; 
Wyatt et al. 1990; De Blasi et al. 1993; Elwell et al. 
1994), To carry out a study, near-infrared light is 
transported to the head by optical fibers which are 
usually terminated with reflective prisms (henceforth 
called optodes) which direct the light into the tissues. 
Light is collected by a similar fiber arrangement and 
detected by a photomultiplier tube or photodiode. The 
separation between input and detection optodes is 
usually in the range 2-4 cm. Changes in optical attenu- 
ation are then measured at a number of wavelengths, 
and these can be related to changes in hemoglobin 
oxygenation and concentration using the known absorp- 
tion spectra of the hemoglobin (Wray et al, 1988). 
Measurements can also be made of time of flight of the 
light or the phase shift of intensity modulated light. A 
review of the different measui einenl types used in NIR 
spectroscopy can be found in Delpy and Cope (1997). 

Recently, several authors (Hoshi and Tamura. 1993; 
Villringer et al.. 1993; Maki et al. 1995; Meek et al. 
1995; Gratton etal. 1995a) reported studies using NIR 
spectroscopy to investigate the response of the brain to 
visual or motor stimuli. For these studies, the optodes 
are placed on the scalp over the appropriate cortical 
area, a suitable stimulus protocol is followed, and 
changes in hemoglobin oxygenation and concentration 
are followed. One study (Gratton et al.. 1995b) also 
reported light-scattering changes in the tissues. Usu- 
ally, the stimulus is repeated several times with a rest 
period between, and, if necessary, the measurements 
are averaged over the stimulus cycle. 

At present, the interpretation of the measured attenu- 
ation is limited by incomplete knowledge of which 
regions in the brain are being sampled by the NIR light, 
although some attempts to define the region have been 
made by comparing NIRS data with that obtained from 
PET (Hoshi et al. 1994; Hock et al.. 1997) and fMRI 
studies (Kleinschmidt et al, 1996). Developments of 
mathematical models of light propagation are ongoing, 
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and recent theoretical modeling (Okada et a]., 1997; 
Firbank et al, 1996; Hielscher et al, 1997) has indi- 
cated that the path taken by the light is greatly affected 
by the CSF, particularly in the adult head. The CSF is a 
clear fluid which cannot be modeled by simple diffusion 
theory approximations of light transport. Results from 
the more recent modeling, which can correctly incorpo- 
rate the effects of the CSF, indicate that the light 
penetration in the adult brain may be limited to the 
outer cortical gray matter (Firbank etal., 1996; Okada 
etal. 1996, 1997). 

The purpose of this investigation, therefore, was to 
model the distribution of light in the brain in a visual 
evoked-response study and to predict what proportion 
of the detected signal arises from cerebral tissue in 
order to help in the interpretation of experimentally 
measured changes in optical attenuation. 

THEORETICAL BASIS OF OPTICAL 
SPECTROSCOPY BV VIVO 

In a clear solution, light attenuation is linearly 
related to the compound concentration, its specific 
absorption coefficient, and the distance traveled through 
the solution. This is known as Beer's law. For a clear 
solution with absorption coefficient Pa in a cuvette of 
thickness d, the attenuation A is given by 

/4 = In (V7) = p,rf, (1) 

where Iq is the incident intensity and / the measured 
intensity. For a scattering medium, the distance the 
light travels is greater than the physical distance 
between the source and the detector, but it has been 
shown (Delpy et al. 1988) that Eq. (1) can still be used 
to relate changes in absorption (Ap^) to changes in 
attenuation {AA), provided that the mean distance, L, 
traveled by the light is known, 

= In = ApaL, (2) 

where Ii and are the intensities measured, and the 
mean distance 

L = ^d, (3) 

P being a numerical factor representing the increase in 
path traveled by the light due to multiple scattering. 

Obviously, Eq. (2) can be rewritten in terms of the 
optical pathlength: 




(4) 



This equation can be used to calculate the optical 
pathlength L in a homogeneous medium from the 
change in attenuation due to a small change in absorp- 
tion coefficient. In the case of a heterogeneous medium, 
where the absorption can change differently in different 
regions, the partial path length can be defined for a 
change in region x. 



Since experimentally it is not possible to vary the 
absorption coefficient of tissues in a controlled manner, 
it is not possible to measure the partial path length of 
different regions of the head in vivo. Measurements in 
phantoms are possible, however, as are calculations of 
the partial path length from accurate mathematical 
models of light transport. 

One implication of the above equation is that if a 
change in absorption (i.e., chromophore concentration) 
occurs in one region only, the change in Pa (and hence 
concentration) can only be calculated if the partial path 
length in that region is known. 

Photon Measurement Density Functions (PMDFs) 

The PMDF describes how sensitive a measurement is 
to different points in the tissue. Two types of measure- 
ment are considered in this paper: simple attenuation 
changes, from which concentration changes can be 
derived as described above, and mean time measure- 
ments. The latter is a measurement of the average time 
taken by light to travel from the source to the detector 
optode. Obviously, this is equivalent to the mean dis- 
tance traveled divided by the speed of light. Mean time 
can be measured either by timing ultrashort pulses of 
light as they pass through the tissue or by using 
amplitude-modulated light, in which case the phase 
shift of the detected light is measured. In tissue, mean 
time (f) and phase shift 4> (in radians) are simply related 
by 




(6) 



provided that the modulation frequency, f, is less than 
200 MHz (Arridge etal. 1992). 

The mean time (or phase shift) can be used as an in 
vivo measurement of total path length. Changes in pa 
and scattering coefficient (ps) will cause changes in it, 
as well as changing attenuation. Gratton et al (1995b) 
reported fast changes (100-200 ms) in the phase shift 
during an evoked-response study which they claim are 
caused by changes in scattering as the neurons depolar- 
ize following the stimulation. 
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The calculation of PMDFs, which has usually been 
done with diffusion theory, is discussed elsewhere 
(Arridge, 1995; Arridge and Schweiger, 1995), but a 
brief definition is as follows. 

For intensity measurements, the PMDF J{p, m, s) at 
position p in tJie tissue, for a source at s and measure- 
ment of change in fluence Ar(m) at position m, is 



Jr{p, m, s) = 



Ar(m) 



Aiia(p) being the change in absorption coefficient at p. 

Similarly, for mean time (f), the corresponding 
PMDF is 



J(4(p, m, s) = 



mm) 



Note that the mean time PMDF can be either positive 
or negative, since an increase in absorption coefficient 
can cause either an increase or decrease in mean time 
depending upon the position of the change relative to 
the source and detector optodes. 

Differential Mean Time Factor 

To calculate the effect of absorption changes in a 
region on the mean time, we can define the differential 
mean time factor (DMTF) for region a: in a similar 
fashion to Eq. (5); 



DTMFW = 



A(f) 



TECHNIQUES 

A MRI scan of one of the authors' heads was taken 
and an axial slice through the region of the visual 
cortex selected. This image was used to generate out- 
lines (Fig. 1) of the gray matter, white matter, CSF, and 
extra cerebral tissue (skull, skin) . These outlines were 
used to generate a finite element mesh (see Fig. 2). The 
approximate radii of the tissue types were also mea- 
sured and used in a separate cylindrical model of the 
tissues. The dimensions and tissue optical properties 
(absorption coefficient p^ and transport scattering coef- 
ficient p J are shown in Table 1. 

The optical properties aie typical of those quoted in 
the literature (Firbank et aL, 1993; van der Zee et al, 
1993; Cheong, 1995) at 800 nm. Oxygenated and deoxy- 
genated hemoglobin have an isobestic point at 798 nm, 
with an absorption coefRcient of 0.2 mm ' mM (Matcher 
et al, 1995). All absorption and scattering coefficient 
are in base e. 

Since diffusion theory is only valid in scattering 
regions, a hybrid diffusion/radiosity model was used to 
calculate light transport in the head. In this scheme, 
which has been specifically developed to cope with 
tissues containing clear regions (Firbank etal, 1996), a 
2D finite element model (FEM) solution of the diffusion 
equation (Arridge etal, 1993) is used to calculate light 
transport in the scattering regions, while in the nonscat- 
tering regions of the CSF (the blank section in the 
middle of the mesh in Fig. 2), radiosity light tracing is 
employed. This hybrid scheme has been validated 
against both experimental measurements in phantoms 
and a Monte Carlo model (Firbank etal., 1996). 

In addition, we performed a Monte Carlo simulation 
(Hlraoka et al., 1993) of a concentric cylinder with 
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FIG. 2. Finite element mesh generated from Fig. 



properties as shown in Table 1 , using 5 million input 
photons. This was a 3D simulation. FEM/radiosity 
calculations were also performed for a 2D cylindrical 
geometry. 

The simplified cylindrical geometry was used in 
order to investigate the effect of changing the thickness 
of the CSF, since it is easier to define a thickness change 
for a uniform layer (in a previous study, we have shown 
that the presence of the sulci have little effect upon the 
thickness of the brain tissue contributing to the NIRS 
signal) (Okada et al. 1997). The FEM model based 
upon the MRI data, however, allowed qualitative analy- 
sis of the effect of the nonuniformity of the gray matter 
surface. The 3D Monte Carlo cylinder model was run to 
investigate the spread of light perpendicular to the 
measurement plane, since this is not possible to do in a 
2D model (which our finite element model is currently 
restricted to because of the computation time required). 

PMDFs were generated with source/detector separa- 
tions of 2, 3, and 4 cm which are typical of the values 
used in the published evoked-response studies. The 
input position for the light in the MRI generated 
geometry is shown in Fig. 2. This was chosen to be 
approximately over the visual cortex while avoiding the 
sagittal sinus. Data for the surface fluence and mean 
time were calculated at all detector positions. 



TABLE 1 





Radius (mm) 


(la (mm"') 


V.s(mm"i) 


White matter 


68.5 


0.005 


6.0 


Gray matter 


75.5 


0.025 


2.5 


CSF 


79.0 


0.003 


0.001 


Skin/skull 


92.0 


0.01 


2.0 



RESULTS 

Figure 3 shows the PMDFs for measurements of 
mean time and intensity for the FEM calculations on 
the head. These show that for a measurement of 
intensity, the region of the brain which Is interrogated 
is limited to the cortical gray matter immediately below 
and between the optodes and that this region is sampled 
approximately uniformly. For a measurement of mean 
time (or phase if using an intensity modulated spectrom- 
eter) a slightly greater region is sampled in terms of 
both depth and lateral spread. The depth of gray matter 
sampled does not increase substantially with increased 
optode spacing. Figures 4a and 4b show profiles of the 
PMDFs in the brain tissue halfway between the opto- 
des and perpendicular to them (in the 2D plane) for 
source/detector spacings of 2, 3, and 4 cm. The profiles 
have been normalized by dividing each of them by the 
received signal at the appropriate spacing. The contri- 
bution of the tissue to the detected signal falls exponen- 
tially with distance into the gray matter; tJie lie 
penetration depth for the intensity PMDF In the gray 
matter can be calculated as 1.5 mm and for the mean 
time as 2.4 mm. 

Figure 5 a shows the optical path lengths in each of 
the tissue types as a function of optode spacing. This is 
proportional to the attenuation change caused by a 
variation in absorption within that region. A 1 pM 
increase in blood concentration in the skin will hence 
give rise to a change in attenuation approximately 20 
times bigger than a 1 \xM change in the gray matter. 
The partial path length through the white matter is 
less than 0.5 mm and is effectively zero on this scale, 
and hence changes in blood concentrations solely occur- 
ring there will not be seen. 
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FIG. 3. Intensity and mean time PMDFs for 2-, 3-, and 4-cm optode separation. For the intensity PMDF, dark color indicates an increase in 
attenuaUon wltli increase in jtg. For the mean time PMDF, dark color Indicates a decrease in mean time with increase in and light color an 
increase in mean time. 



Figure 5b shows a similar graph for the mean time, 
the y axis showing the change in mean time for a 
change in absorption coefficient. The mean time mea- 
surement is approximately twice as sensitive to changes 
occurring in the brain than the intensity measurement, 
although it still is a factor of 10 less than the surface 
regions. 

Information about the spread of photons perpendicu- 
lar to the source/detector plane was gained from the 3D 
concentric cylinder Monte Carlo model by recording the 
position of all detected photons as they passed through 
the surface of the gray matter. The resulting intensity 



PMDFs are shown in Fig. 6. Figure 7 shows the profile 
through the PMDFs along the axis indicated by the 
dotted line normalized by the total detected intensity. 
As can be seen, the Mght probes an area of brain 
approximately 1 cm on either side of the source/ 
detector plane. This vertical spread does not change 
significantly with source/detector spacing. It is worth 
noting that the use of a cylindrical as opposed to a 
spherical 3D model would, if anything, slightly overes- 
timate the vertical spread to be expected, so the values 
presented here probably represent a maximum extent 
to the sampled tissue volume (Okadaefa7., 1996, 1997). 
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FIG. 4. Profile through (a) intensity and (b) mean time PMDFs in gray n- 
calculation on the MRI-generated mesh. 
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Effect of Brain Movement 

In NIRS studies, optical changes can be observed in 
synchrony with both the heart and the respiratory 
rates, and these need to be filtered out in measuring the 
small evoked-response signal (Gratton and Corballis, 
1995). These signals have usually been assumed to be 
due to changes in blood flow/volume that occur with 
each heart beat or respiration (Elwell etal, 1994), but 
it is also possible that they can arise due to movements 
of the brain surface or from pial cerebral vessels lying 
on the brain surface. The brain is not fixed inside the 
skull and can move or expand to a limited degree. The 
pulsatile nature of the blood flow causes the brain to 



throb with the heart beat. This has been observed in 
patients undergoing brain surgery, although in NIR 
studies, of course, the skull is not opened, and the brain 
has less room for movement. 

Changes in the blood volume of the brain (as part of 
an evoked response) will also lead to a change in brain 
volume and will cause it to expand or shrink, displacing 
CSF. Both the absorption and the volume change may 
alter the light distribution in the head, resulting in an 
optically observable signal. We have applied our previ- 
ously described modeling to determine the relative 
contributions of movement and absorption variation to 
the NIRS signal. 




Source/detector separation (mm) 

FIG. 5. (a) Partial path lengths though skin/skull, CSF, and gray matter for FEM/radiosity calculation on the MRI-generated mesh and 
Monte Carlo modeling of the 3D cylinder and (b) differential mean time factor for different regions for FEM/radiosity calculation on the 
MRI-generated mesh. 
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FIG. 7. Profile of detected photons as they cross the gray matter 
boundary for 2-, 3-, and 4 cm source/detector spacing. 



Greitz et al. (1992) reported pulsatile brain move- 
ment using magnetic resonance phase imaging in adult 
volunteers. This movement was greatest around the 
base of the brain (0. 1 mm displacement and velocities of 
up to 2 mm/s) . The magnitude of the motion decreased 
toward the periphery of the brain, however, and was 
less than 0.01 mm at the brain surface. 

Figure 8 shows the calculated change in attenuation 
and mean time for a decrease of 0, 1 mm of the CSF 
thickness in the simplified cylindrical model (i.e., the 
inner radius of the CSF changes from 75.5 to 75.6 mm, 
with all other parameters remaining constant). 

When the blood volume in Lhe brain increases, the 
attenuation will change for two reasons: (a) because the 
blood content has increased, causing higher absorption, 
and (b) because the brain will expand slightly, changing 
the distribution of light in the head. Figure 9 shows a 



Angle(degrees) 

FIG. 6. Monte Carlo-generated intensity PMDFs for the gray 
matter/CSF interface, showing the vertical spread of the tight on the 
gray matter surface. PMDFs are shown for light detected at 2-, 3-, 
and 4 cm separation. The positions of the optodes are shown with a 
cross. The dotted line shows the directloin across which the profiles 
were generated (see Fig. 7). 




20.0 30.0 40.0 
Positiofi (nnni} 

FIG. 8. Change in attenuation and mean time for CSF thickness 
decreasing from 3.5 to 3.4 mm (in 2D cylinder FEM/radioslty). 
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FIG. 9, Comparison of signal change for (a) attenuation and (b) mean time caused by an absorption increase with that caused by brain 
expansion due to increased blood volume (from 2D cylinder FEIVl/radiosity calculation). 



comparison of these two effects, i.e., an increased 
absorption, but constant geometry, and an increased 
brain volume at constant absorption. The data are 
calculated for a 1 liM change of hemoglobin concentra- 
tion in the gray matter. The gray matter thickness is 7 
mm and the brain blood volume is approximately 5% 
with a typical hemoglobin concentration of 84 jiM 
(Sakai et al, 1985), so a 1 pM increase is equivalent to 
5 X 1/84 = 0,06% increase in volume and so assuming 
that the brain swells outward only and that the brain 
hematocrit stays constant, this is equivalent to an 
increase of 0.004 mm in thickness. 

IMPLICATIONS 

When used on the scalp surface, NIR spectroscopy 
can pick up changes in hemoglobin concentration on 
the cortical surface. The signal is probably limited to 
changes in hemodynamics occurring in the top 2-3 mm 
of the cortex. The area of the brain surface which 
contributes to the signal is approximately as wide as 
the optode spacing and extends laterally 1 cm either 
side perpendicular to the optode position. In NIR 
spectroscopy, the calculated quantity ACL is a product 
of the changes in concentration and optical pathlength, 
To derive absolute changes in concentration, the optical 
path length must be known. If the concentration changes 
are limited to one region only, the partial path length in 
that region must be known to calculate the change in 
concentration correctly. In the case of evoked response 
studies, the hemoglobin concentration changes are 
assumed to occur in the brain only. Although the 
PMDFs are weighted toward surface tissues, this does 
not necessarily mean that NIRS signal comes mainly 
from the surface. This is because in evoked-response 
studies, we are measuring changes in signal, and if the 



only change in hemodynamics occurs in the brain, then 
that is what will be detected. The technique is, how- 
ever, more sensitive to changes occurring in the skin 
and skull than in the brain. 

While the exact dependence of the NIR signal contri- 
bution in an individual study will depend upon the 
geometry and the exact optical properties, in our theo- 
retical model, we found that a 1 pM change in blood 
concentration will give rise to approximately 0.001 OD 
change in attenuation at 4-cm spacing and approxi- 
mately 1-ps change in mean time. Typical instrumental 
noise at 3-cm spacing in tissue is 0.001 OD for attenua- 
tion and 0. 1 " in phase (equal to 1 . 4 ps at 200 MHz) . This 
tji^ically gets a factor of 3 worse for every centimeter 
increase in optode spacing and is essentially limited by 
the detected intensity and the quantum efficiency of the 
detecting device. The signal can of course be improved 
by averaging and performing several studies on the 
same volunteer. 

It is interesting to note that, although experimen- 
tally the NIRS signal changes measured over the visual 
cortex vary considerably in magnitude from subject to 
subject, the average level of intensity change (Meek et 
al, 1995) and phase change (Gratton et a!.. 1995a) 
observed is approximately five times greater than that 
predicted from this simple modeling in which the 
hemoglobin changes are assumed to occur only in the 
brain parenchyma. We have previously speculated 
(Okada et al, 1997, 1998) that the considerable time 
spent by the light in the clear CSF would lead to it 
contributing significantly to the NIRS signal were an 
absorption change to occur within it. In reality, the CSF 
layer comprises the subarachnoid space which includes 
pial cerebral vessels lying on the brain surface. Any 
small change in the hemoglobin concentration, oxygen- 
ation, or size of these vessels as a result of an evoked 
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response is likely therefore to be a major contributor to 
the experimentally obsei'ved signal. The current hybrid 
diffusion/radiosity model cannot deal with this geom- 
etry, but work is under way to modify it to enable the 
effect of discrete absorbers within the clear layer to be 
modeled. 
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ABSTRACT 

The results of measurements of transmittance of higfi power laser irradiation throu^ skull bones and scalp are 
presented. Dependences of transmittance on sample thickness were received. Character of transmittance was 
investigated and characteristics of heterogeneity of the scattering stnicture of the skull bones are proposed. Besides that, 
variation of temperature of skull and scalp surfaces under exposure of high power laser irradiation during experiments 
was controlled. Experimental results were verified by Monte-Carlo simulations. 

Keywords: transmittance, bone structure, laser therapy, Monte-Carlo simulation 

1. INTRODUCTION 

During the last years laser technique had found increasingly wide application in medicine. In addition to using of hi^- 
power lasers in surgery (laser scalpel) laser irradiation is widely used in therapy. There are a lot of fields of application 
of laser therapy, such as: cure of arthritb, oculars, cardiac, skins and blood diseases. The influence of different types of 
laser systems on different biological otgects was investigated during this time. 

LasCT therapy at present is one of the most effective methods of treatment and prophylaxis of various diseases, primarily 
owing to possibility of fine adjustment of biological processes and noninvasive influence on structure of biological 
tissues. However, it is necessary to define properly radiation dose to avoid of negative thermal effect on biological 
tissues. As is well known, the power density of laser irradiation within the range of from 0. 1 to 10 mW/cm2 is had the 
most effective influence to human tissues. The lower power density of laser irradiation could exert just stimulation 
therapy; the abnormal changes in tissues could h^yipen using the high value of power dcaisity of laser irradiation. 

One of the possible adaptations of this method is a therapy of stroke and in&rcts. Accordingly, the main aim of our 
studies was estunation of a dose of laser irradiation transmitted through skull bones and scalp tissue and arrived at a 
brain. A lot of parameters are significant in this case such as wavelength and power of iiradiation, geometry of a beam, 
aind mode of irradiation (continues wave, pulsed beam). To understand the adaptation processes within bone dssue it is 
not only necessary to take into consideration bone density but also bone stiucUire. 

2. EXPERIMENTS 

Experimental research was carried out with preparing bone tissue and skin of the head. The scheme of the experimental 
setup is presented in Figure la. Coltim^ed, S 10 nm, 1 Watt laser beam of approximately 30 mm in diameter was used 
as a light source. 

CCD cam^a with necessary neutral filters was used for obtaining of a pattern of spafial distribution of transmitted 
irradiation. 
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Figure I; Experimental setup and method of measurements 



The absolute value of power at the central of distribution (maximum power) was tneasured using a 1 mm diameter fiber 
with cosiiie corrector and an Ocean Optics a USB2000 spectrophotometer. Thickness of the sample was measured at 
five points witli the help of the shde gauge and averaged. The skin samples were held between two object-plates during 
measurement. Temperature of samples was measured by the voltmeter V7-27A/1 with tiie contact probe for temperature 
detection (Fig. lb). Square of working substance of the probe is equal to 1 square nuUfmeter. Temperature was 
measured at the skin surfece and skull in the center of uradiated ft eld after exposure. 

3. MEASUREMENTS 



Using CCD we made a sequence of frames (25 frames) with each domain of interest (temples, forehead, occiput, and 
vertex) for averaging, sequence of frames with the laser switch off for subtraction of backgrounds, and one frame wifli 
the scale for calculation of real dimensions of the frames. Then maximum power in tiie center of the distribution was 
measured five times by the spectrometer, using available calibration, and the light distribution was normalized to mean 
value. TTic resulting power density distribution was calculated in ttie following way: 

( Pdistr^-\- Pdistr--{r Pdistr2^ Pbackgnd^ + ...+ Pbackgnd2^\ (i'maxl + ...+ PmaxS) 
"'^'^^ ' { 25 25 J SPdistrmsx 

Hie same procedure as described above was realized both with the skull only and witii the skull with shaved scalp, Afier 
measur^ents we calculated the main characteristics of light transmittance: average radius of light distribution with tiie 
fixed power density, and its standard deviation. For example: 



^0.7max ^ ' 

where is a point from vicinity of 0,7MaxPower level (0.7-e; 0,7+e), and /?^„ is a radius vector of this point on flie 
)tion of ^0 is a cent» of coordinate. amount of this points. Analogously for Ro.s and R^t. 



4. RESULTS AND CONCLUSIONS 
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Measurements were carried out on 20 cadaver's heads, at 5 domains: foreliead, left and right temples, occ^)ut and 
vertex, bi such a way dependence of transmittance on sample thickness was build up on 100 points. Measurements were 
performed both on the skull with shaved scalp and on the separate skull. In Figure 2 dqranieaice of transmittance on 
sample thickness for bo& cases is presented. 

As it can be seen fi'om this picture, decrease of transmittance with increase of sample thickness can be approximated 
well by expon^ial decay. The value of transmittance is varied from 0.5% to 5% in case of sample with scalp and from 
1% to 16% in case of single skull. In consequence of input power density was 68 mW/cm^ value of transmitted power 
deosify varied from 340 (iW/on^ to 3400 (iW/ctn^ and from 680 )iW/cm^ to 9000 fiW/cm' respectively. It is eiddent 
from results of experiments with separated bone tissue and experiments with skull and scalp, that layer of skia 
influences on transmittance the same way as skull, with that thicloiess of skin layer is less tiian thickness of skufL Such 
experimental results comply well with theoretical Monte-Carlo simulations, results of which are presented in Figure 3. It 
is obvious from this picture, that the theoretical dependence ol transmiltance on sample thickness has the same decay 
rate as exp^imental curve. Distinctions of the theoretical values of transmittance of skull-scalp sample from 
experimental data can be explained by that average optical parameters were used in Monte-Carlo simulations. 

Temperature control was provided during experiments. Temperature of the skin surface rose on 4 °C on average for the 
time of exposure of 4 minutes and the temperature of the externa! surface of the skull rose on 2 °C on average for the 
same time of exposure under ^uch power of irradiation. Negligible changes of the temperahire of the internal sur&ce of 
the skull in few percent of degree were notice. It can be referred to inaccuracy of measurranents. Dependence of 
tenq)erature on time of exposure for one of the samples is represented in Table 1. 




Figure 2: Dependence of transmittance on thickness of the sample - experimental resufts 

It was said above, that the radii of equal power density (0.1, 0.5 and 0.7 of maximum value) were estimated for the 
distributions of power density obtained in experiments. The dependences of radius 0.5 of maumum on transmittance are 
presented in Figure 4. It follows from this ^gure, that increasing of radius O.Smaximum with increasing of transmittance 
is well seen in most cases. Besides that, standard deviation of each radius of fixed power density was calculated, which 
defined, in this case, degree of heterogeneity of tissue's transmittance related witii the features of the tissue structure, 
maii^y with the features of the bone tissue. 

We can suppose that this scheme of irradiation can be used for therapeutic irradiation and can be employed for clinical 
study on the assumption of expected therapeutic power density at the dura maier of 7.5 mW/cm^. Values of power 
density obtained in both experiment and simulation lie in titis limit. 
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Figure 3: Dependence of transmittance on thickness of sample - Monte^Carlo simulations 



Table I: Dependence of temperature on time of exposure 


object 


time of irradiation, min 


temperature, °C 


skull 


1 


21.9 


scalp 


1 


22.0 


skull 


2 


22.6 


scalp 


2 


23.2 


skull 


3 


24.1 


scalp 


3 


24.6 


Skull 


4 


24.8 


scalp 


4 


25.6 




I"*™*- * tiansmittanca, % 

Figure 4: Depraidences of radius 0.5 of maximum on transmittance 
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